Conteúdos abordados neste laboratório

• Primeiros passos no R e no Rstudio
• Explorando uma base de dados carregada da nuvem
• fazendo filtros e gráficos com o R

Carregamento de pacotes

O R é uma linguagem de programação de código aberto, o que significa que os usuários podem contribuir com pacotes que facilitam nossas vidas. Ainda melhor, podemos usá-las gratuitamente. Para este laboratório, e muitos outros futuros, usaremos principalmente os seguintes pacotes do R:

• dplyr: para descrição de dados
• ggplot2: para visualização de dados
• kableExtra: Para deixar as saídas de dados em formato de tabela
• rmarkdown: para fazer relatórios como este, usando markdown

Para carregá-los no nosso ambiente de trabalho, digitamos no console:

library(ggplot2)
library(plyr)
library(tidyr)
library(dplyr)
library(DescTools)
library(htmltools)
library(knitr)
library(rmarkdown)
library(combinat)
library(kableExtra) #Para gerar tabelas com um layout agradável

Os dados: Registros de batismos do Dr. Arbuthnot

Aqui, vamos carregar os dados diretamente da itnernet, na página do OpenIntro Statistics:

source("http://www.openintro.org/stat/data/arbuthnot.R")

O conjunto de dados Arbuthnot se refere ao Dr. John Arbuthnot, um médico, escritor e matemático do século 18. Ele se interessou pela razão ente o número de meninos e meninas recém-nascidos, e para isso ele coletou os registros de batismo de crianças nascidas em Londres no período compreendido entre 1629 e 1710. Vamos começar dando uma olhada nos dados, com o auxílio da função head():

arbuthnot %>%
  head(10) %>%
  kable() %>%
  kable_styling(
    full_width = F, 
    bootstrap_options = c("striped", "hover", "condensed",
                          "responsive")
    )

Podemos ver as dimensões desse banco de dados

dim(arbuthnot)

## [1] 82  3

Para ver os nomes das variáveis observadas (colunas), fazemos:

names(arbuthnot)

## [1] "year"  "boys"  "girls"

Explorando os dados

Nós podemos acessar separadamente os dados de uma única coluna da base de dados da seguinte forma:

arbuthnot$boys

##  [1] 5218 4858 4422 4994 5158 5035 5106 4917 4703 5359 5366 5518 5470 5460 4793
## [16] 4107 4047 3768 3796 3363 3079 2890 3231 3220 3196 3441 3655 3668 3396 3157
## [31] 3209 3724 4748 5216 5411 6041 5114 4678 5616 6073 6506 6278 6449 6443 6073
## [46] 6113 6058 6552 6423 6568 6247 6548 6822 6909 7577 7575 7484 7575 7737 7487
## [61] 7604 7909 7662 7602 7676 6985 7263 7632 8062 8426 7911 7578 8102 8031 7765
## [76] 6113 8366 7952 8379 8239 7840 7640

arbuthnot$girls

##  [1] 4683 4457 4102 4590 4839 4820 4928 4605 4457 4952 4784 5332 5200 4910 4617
## [16] 3997 3919 3395 3536 3181 2746 2722 2840 2908 2959 3179 3349 3382 3289 3013
## [31] 2781 3247 4107 4803 4881 5681 4858 4319 5322 5560 5829 5719 6061 6120 5822
## [46] 5738 5717 5847 6203 6033 6041 6299 6533 6744 7158 7127 7246 7119 7214 7101
## [61] 7167 7302 7392 7316 7483 6647 6713 7229 7767 7626 7452 7061 7514 7656 7683
## [76] 5738 7779 7417 7687 7623 7380 7288

Nesse caso, os dados são impressos em tela como vetores, e não como um banco de dados estruturado.

O R tem algumas funções poderosas para criar gráficos. Podemos criar um gráfico simples do número de meninas batizadas por ano com o comando

arbuthnot %>%
  ggplot(aes(x = year, y = girls)) + 
  geom_point(size = 2, col = "red") +
  labs(
    title = "Registros de nascimento do dr. Arbuthnot"
  ) +
  theme(text = element_text(size = 16))

Podemos também fazer um gráfico de linhas, ao invés de plotar os pontos dispersos:

arbuthnot %>%
  ggplot(aes(x = year, y = girls)) + 
  geom_line(linewidth = 1, col = "red") +
  labs(
    title = "Registros de nascimento do dr. Arbuthnot"
  ) +
  theme(text = element_text(size = 16))

Há alguma tendência aparente no número de meninas batizadas ao longo dos anos? Como você a descreveria?

Para saber o número total de crianças batizadas em cada ano, fazemos:

arbuthnot$boys + arbuthnot$girls

##  [1]  9901  9315  8524  9584  9997  9855 10034  9522  9160 10311 10150 10850
## [13] 10670 10370  9410  8104  7966  7163  7332  6544  5825  5612  6071  6128
## [25]  6155  6620  7004  7050  6685  6170  5990  6971  8855 10019 10292 11722
## [37]  9972  8997 10938 11633 12335 11997 12510 12563 11895 11851 11775 12399
## [49] 12626 12601 12288 12847 13355 13653 14735 14702 14730 14694 14951 14588
## [61] 14771 15211 15054 14918 15159 13632 13976 14861 15829 16052 15363 14639
## [73] 15616 15687 15448 11851 16145 15369 16066 15862 15220 14928

Nesse caso, visualizamos um vetor com 82 entradas, cada uma contendo a soma do número de meninos e de meninas batizados no respectivo ano. Nós podemos ainda utilizar este vetor para gerar alguns gráficos, então iremos salvá-lo como uma coluna permanente do nosso banco de dados:

arbuthnot = mutate(arbuthnot, total = boys + girls)

arbuthnot %>%
  head(10) %>%
  kable() %>%
  kable_styling(
    full_width = F, 
    bootstrap_options = c("striped", "hover", "condensed",
                          "responsive")
    )

Assim, uma nova coluna chamada “total” foi incluída no banco de dados. Podemos criar um gráfico do número total de batismos por ano com o comando:

arbuthnot %>%
  ggplot(aes(x = year, y = total)) + 
  geom_line(linewidth = 1, col = "orange") +
  labs(
    title = "Registros de nascimento do dr. Arbuthnot"
  ) +
  theme(text = element_text(size = 16))

Podemos também calcular a proporção do número de meninos batizados em relação ao número de meninas batizadas em cada ano:

arbuthnot$boys / arbuthnot$girls

##  [1] 1.114243 1.089971 1.078011 1.088017 1.065923 1.044606 1.036120 1.067752
##  [9] 1.055194 1.082189 1.121656 1.034884 1.051923 1.112016 1.038120 1.027521
## [17] 1.032661 1.109867 1.073529 1.057215 1.121267 1.061719 1.137676 1.107290
## [25] 1.080095 1.082416 1.091371 1.084565 1.032533 1.047793 1.153901 1.146905
## [33] 1.156075 1.085988 1.108584 1.063369 1.052697 1.083121 1.055242 1.092266
## [41] 1.116143 1.097744 1.064016 1.052778 1.043112 1.065354 1.059647 1.120575
## [49] 1.035467 1.088679 1.034100 1.039530 1.044237 1.024466 1.058536 1.062860
## [57] 1.032846 1.064054 1.072498 1.054359 1.060974 1.083128 1.036526 1.039092
## [65] 1.025792 1.050850 1.081931 1.055748 1.037981 1.104904 1.061594 1.073219
## [73] 1.078254 1.048981 1.010673 1.065354 1.075460 1.072132 1.090022 1.080808
## [81] 1.062331 1.048299

Podemos também calcular a proporção de meninos batizados ano a ano em relação ao total de crianças batizadas no respectivo ano:

arbuthnot$boys / (arbuthnot$boys + arbuthnot$girls)

##  [1] 0.5270175 0.5215244 0.5187705 0.5210768 0.5159548 0.5109082 0.5088698
##  [8] 0.5163831 0.5134279 0.5197362 0.5286700 0.5085714 0.5126523 0.5265188
## [15] 0.5093518 0.5067868 0.5080341 0.5260366 0.5177305 0.5139059 0.5285837
## [22] 0.5149679 0.5322023 0.5254569 0.5192526 0.5197885 0.5218447 0.5202837
## [29] 0.5080030 0.5116694 0.5357262 0.5342132 0.5361942 0.5206108 0.5257482
## [36] 0.5153557 0.5128359 0.5199511 0.5134394 0.5220493 0.5274422 0.5232975
## [43] 0.5155076 0.5128552 0.5105507 0.5158214 0.5144798 0.5284297 0.5087122
## [50] 0.5212285 0.5083822 0.5096910 0.5108199 0.5060426 0.5142178 0.5152360
## [57] 0.5080788 0.5155165 0.5174905 0.5132301 0.5147925 0.5199527 0.5089677
## [64] 0.5095857 0.5063659 0.5123973 0.5196766 0.5135590 0.5093183 0.5249190
## [71] 0.5149385 0.5176583 0.5188268 0.5119526 0.5026541 0.5158214 0.5181790
## [78] 0.5174052 0.5215362 0.5194175 0.5151117 0.5117899

Vamos usar novamente a função mutatepara incluir uma nova coluna em nosso conjunto de dados para a proporção do número de meninos em relação ao total de recém-nascidos em cada ano:

arbuthnot = mutate(arbuthnot, ratio_boys = arbuthnot$boys / arbuthnot$total)

arbuthnot %>%
  head(10) %>%
  kable() %>%
  kable_styling(
    full_width = F, 
    bootstrap_options = c("striped", "hover", "condensed",
                          "responsive")
    )

Podemos também fazer um gráfico da proporção de meninos batizados relativo ao total ano a ano:

arbuthnot %>%
  ggplot(aes(x = year, y = ratio_boys)) + 
  geom_line(linewidth = 1, col = "skyblue") +
  labs(
    title = "Porporção de meninos batizados"
  ) +
  theme(text = element_text(size = 16))

Nesse caso, percebe-se que a proporção de meninos oscila em torno de 50%.

Podemos ainda criar uma variável lógica que indica se o número de meninos nascidos é maior que o número de meninas:

arbuthnot = mutate(arbuthnot, more_boys = boys > girls)

arbuthnot %>%
  head(10) %>%
  kable() %>%
  kable_styling(
    full_width = F, 
    bootstrap_options = c("striped", "hover", "condensed",
                          "responsive")
    )

Um novo conjunto de dados: registros atuais de nascimento nos Estados Unidos

Vamos agora trabalhar com um conjunto de dados obtido mais recentemente. Primeiro, carregamos o dataframe:

source("http://www.openintro.org/stat/data/present.R")

1. Quais anos estão incluídos nesse conjuntos de dados?

names(present)

## [1] "year"  "boys"  "girls"

present$year

##  [1] 1940 1941 1942 1943 1944 1945 1946 1947 1948 1949 1950 1951 1952 1953 1954
## [16] 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969
## [31] 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984
## [46] 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999
## [61] 2000 2001 2002

Vemos então que os anos incluídos vão de 1940 a 2002.

2. Quais as dimensões da base de dados?

dim(present)

## [1] 63  3

3. Quais sãos os nomes das colunas ou variáveis?

names(present)

## [1] "year"  "boys"  "girls"

4. Como esses dados se comparam aos de Arbuthnot? Eles estão em uma escala similar?

present %>%
  ggplot(aes(x = year, y = girls)) + 
  geom_point(size = 2, col = "red") +
  labs(
    title = "Registros de nascimento atuais"
  ) +
  theme(text = element_text(size = 16))

cat("Intervalo de anos para os dados do dr. Arbuthnot:", 
    c(min(arbuthnot$year), max(arbuthnot$year)),
    "\nIntervalo de anos para os dados atuais:",
    c(min(present$year), max(present$year))
    )

## Intervalo de anos para os dados do dr. Arbuthnot: 1629 1710 
## Intervalo de anos para os dados atuais: 1940 2002

Como podemos ver a partir do resultado do gráfico acima, bem como com a saída da função cat(), os dados dos conjuntos arbuthnot e present estão em escalas diferentes. No caso dos dados do dr. Arbuthnot, o número de batismos de menina, por exemplo, está na escala de alguns milhares, ao passo que no conjunto de dados atual, esse número está na escala de milhão. Além disso, o intervalo temporal de coleta dos dois conjuntos de dados são também diferentes.

5. Crie um gráfico que mostre a proporção de meninos em relação às meninas para cada ano do conjunto de dados. O que você vê? A observação de Arbuthnot de que os meninos nascem numa proporção maior do que as meninas se mantém nos EUA?

present = mutate(present, ratio_boys_girls = present$boys/present$girls)

present %>%
  ggplot(aes(x = year, y = ratio_boys_girls)) + 
  geom_line(linewidth = 1, col = "orange") +
  labs(
    title = "Razão entre o número de meninos e de meninas",
    subtitle = "Conjunto de dados atual."
  ) +
  theme(text = element_text(size = 16))

Conforme podemos ver na saída do comando abaixo, a observação de Arbuthnot se mantêm nos EUA, embora o gráfico acima nos mostre que a razão entre o número de meninos e de meninas nascidos vem caindo ao longo das décadas, se aproximando de 1.

present$meninos > present$meninas

## logical(0)

Outra maneira de ver se existe algum ano no qual o número de meninos nascidos é menor do que o número de meninas é fazendo um filtro no dataframe:

present[present$boys < present$girls, ]

## [1] year             boys             girls            ratio_boys_girls
## <0 linhas> (ou row.names de comprimento 0)

Como vemos, a saída é vazia, tendo apenas os nomes das variáveis, sem qualquer observação nas linhas.

6. Em qual ano se verifica o maior número de nascimentos nos EUA? Dica: primeiramente calcule o total e salve em uma nova variável. Então, ordene seu banco de dados em ordem decrescente na colula total. Você pode fazer isso interativamente no seu Data Viewer, clicando nas flechas próximas aos nomes das variáveis.

Primeiro, vamos acrescentar a coluna “Total” ao nosso dataframe.

present = mutate(present, Total = present$boys + present$girls)

Em seguida, vamos aplicar um filtro ao dataframe, assim como fizemos acima, e mostrar todas as colunas do dataframe, porém apenas a linha correspondente ao valor máximo da coluna Total:

present[present$Total == max(present$Total), ] 

##    year    boys   girls ratio_boys_girls   Total
## 22 1961 2186274 2082052         1.050057 4268326

Uma maneira mais elegante de fazer esse filtro é por meio da função subset():

present%>%
  subset(Total == max(Total))

##    year    boys   girls ratio_boys_girls   Total
## 22 1961 2186274 2082052         1.050057 4268326

Por outro lado, para responder à mesma pergunta e mostrar apenas as colunas ano e Total, na linha correspondente ao total máximo, fazemos:

present[present$Total == max(present$Total), c('year', 'Total')]

##    year   Total
## 22 1961 4268326

Ou, por meio da função select():

present%>%
  subset(Total == max(Total), select = c("year", "Total"))

##    year   Total
## 22 1961 4268326

Esses dados são provenientes de uma pesquisa realizada pelo Centro de Controle de Doenças (Center For Disease Control) Link do Artigo.

Conteúdos abordados neste laboratório

• Histogramas
• Estatísticas descritivas
• Tabelas de contingência
• Aprofundando a biblioteca dplyr
• Observando relações entre duas variáveis quantitativas
• Observando relações entre duas variáveis qualitativas
• Observando relações entre uma variável qualitativa e uma quantitativa
• Box plots
• Distribuição normal

Carregamento e detalhamento dos dados

O Sistema de Monitoramento de Fatores de Risco Comportamental (Behavioral Risk Factor Surveillance System, BRFSS) é uma pesquisa anual realizada por telefone com 350.000 pessoas nos Estados Unidos. Como seu nome implica, o BRFSS foi desenvolvido para identificar fatores de risco na população adulta e relatar tendências emergentes na saúde. Por exemplo, os respondentes são indagados sobre sua dieta e atividades físicas semanais, seu diagnóstico de HIV/AIDS, uso provável de tabaco, e mesmo seu nível de cobertura por planos de saúde. O website do BRFSS http://www.cdc.gov/brfss contém uma descrição completa desta pesquisa, incluindo as questões que motivaram o estudo e muitos resultados interessantes derivados dos dados.

Nós focaremos em uma amostra aleatória de 20.000 pessoas do BRFSS conduzida em 2000. Ainda que existam mais de 200 variáveis neste conjunto de dados, nós trabalharemos com um subconjunto menor.

Começamos importando os dados das 20.000 observações para dentro do espaço de trabalho do R. Depois de inicializar o RStudio, entre com o seguinte comando:

source("http://www.openintro.org/stat/data/cdc.R")

O conjunto de dados cdc que aparece em seu espaço de trabalho é uma matriz de dados, onde cada linha representa um caso e cada coluna representa uma variável. Lembre-se que o R denomina este formato de dados como banco de dados (data frame), que será um termo utilizado ao longo dos laboratórios.

Para visualizar o nome das variáveis, digite o comando:

names(cdc)

## [1] "genhlth"  "exerany"  "hlthplan" "smoke100" "height"   "weight"   "wtdesire"
## [8] "age"      "gender"

Ainda, nosso banco de dados tem 20000 indivíduos entrevistados, sendo 10431 do gênero feminino.

Cada uma dessas variáveis corresponde a uma questão que foi feita na pesquisa. Por exemplo, para genhlth, os respondentes foram indagados sobre sua saúde geral, respondendo good, very good, excellent e fair. A variável exerany indica se o respondente se exercitou no último mês (1) ou não (0). Da mesma forma, hlthplan indica se o respondente tem alguma forma de cobertura por plano de saúde (1) ou não (0). A variável smoke100 indica se o respondente fumou pelo menos 100 cigarros ao longo da vida. As variáveis height, weight, wtdesire, age e gender registram a altura (em polegadas), o peso (em libras), o peso desejado (em libras), a idade (em anos) e o gênero, respectivamente.

Em tese, você poderia verificar toda a base de dados de só uma vez digitando seu nome no console, mas isso pode não ser muito sábio neste caso. Sabemos que cdc tem 20.000 linhas, portanto verificar o conjunto de dados inteiro significa inundar sua tela. Uma função bastante útil para dar uma rápida olhada no seu conjunto de dados, ver suas dimensões e os tipos de dados é a função str:

str(cdc)

## 'data.frame':    20000 obs. of  9 variables:
##  $ genhlth : Factor w/ 5 levels "excellent","very good",..: 3 3 3 3 2 2 2 2 3 3 ...
##  $ exerany : num  0 0 1 1 0 1 1 0 0 1 ...
##  $ hlthplan: num  1 1 1 1 1 1 1 1 1 1 ...
##  $ smoke100: num  0 1 1 0 0 0 0 0 1 0 ...
##  $ height  : num  70 64 60 66 61 64 71 67 65 70 ...
##  $ weight  : int  175 125 105 132 150 114 194 170 150 180 ...
##  $ wtdesire: int  175 115 105 124 130 114 185 160 130 170 ...
##  $ age     : int  77 33 49 42 55 55 31 45 27 44 ...
##  $ gender  : Factor w/ 2 levels "m","f": 1 2 2 2 2 2 1 1 2 1 ...

Note que o R chama as as variáveis categóricas de factor e a quantidade de categorias de levels. Você também poderá dar pequenas espiadas nos dados utilizando os comandos head (que mostra as primeiras linhas de dados). Vamos utilizar também a função kable(), do pacote kableExtra, a fim de deixar a tabela com um layout agradável:

head(cdc) %>%
  kable() %>%
  kable_styling(
    full_width = F, 
    bootstrap_options = c("striped", "hover", "condensed",
                          "responsive")
    )

tail() (que mostra as últimas linhas de dados):

tail(cdc) %>%
  kable() %>%
  kable_styling(
    full_width = F, 
    bootstrap_options = c("striped", "hover", "condensed",
                          "responsive")
    )

ou as técnicas de construção de subconjuntos que você aprenderá no decorrer do laboratório.

Por outro lado, se quisermos mostrar o banco de dados na saída do arquivo html como uma tabela estruturada, podemos usar a biblioteca rmarkdown:

paged_table(cdc)

Podemos também estruturar a tabela de saída por meio da biblioteca kableExtra, mudando o nome das colunas da maneira que achar conveniente:

kable(head(cdc, 10), col.names = c("Saúde", "Atividade Física", "Plano de Saúde", "Fumou nos últimos 100 dias?", "Altura", "Peso", "Peso desejado", "Idade", "Gênero")) %>% 
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed", "responsive"))

Descrevendo Distribuições

cdc %>%
  ggplot(aes(x = age)) + 
  geom_histogram(col = "black", fill = "orange") + 
  labs(
    x = "idade",
    y = "contagem",
    title = "Histograma da idade"
  ) +
  theme(text = element_text(size = 18))

cdc %>%
  ggplot(aes(x = age)) + 
  geom_histogram(col = "black", fill = "orange", binwidth = 1) + 
  labs(
    x = "idade",
    y = "contagem",
    title = "Histograma da idade"
  ) +
  theme(text = element_text(size = 18))

cdc %>%
  ggplot(aes(x = age)) + 
  geom_histogram(col = "black", fill = "orange", binwidth = 20) + 
  labs(
    x = "idade",
    y = "contagem",
    title = "Histograma da idade"
  ) +
  theme(text = element_text(size = 18))

cdc %>%
  ggplot(aes(x = age)) + 
  geom_histogram(col = "black", fill = "orange", alpha = .4, bins = 30) + 
  labs(
    x = "idade",
    y = "contagem",
    title = "Histograma da idade"
  ) +
  theme(text = element_text(size = 18))

cdc %>%
  ggplot(aes(x = age)) + 
  geom_histogram(
    aes(y = after_stat(density)), 
    col = "black", 
    fill = "orange", 
    alpha = .4, 
    bins = 30
    ) + 
  labs(
    x = "idade",
    y = "contagem",
    title = "Histograma da idade"
  ) +
  theme(text = element_text(size = 18))

cdc %>%
  ggplot(aes(x = age)) + 
  geom_histogram(
    aes(y = after_stat(density)), 
    col = "black", 
    fill = "orange", 
    alpha = .4, 
    bins = 30
    ) + 
  geom_density(
    col = "red", 
    linewidth = 1
    ) +
  labs(
    x = "idade",
    y = "contagem",
    title = "Histograma da idade"
  ) +
  theme(text = element_text(size = 18))

cdc %>%
  select(weight) %>%
  summary()

##      weight     
##  Min.   : 68.0  
##  1st Qu.:140.0  
##  Median :165.0  
##  Mean   :169.7  
##  3rd Qu.:190.0  
##  Max.   :500.0

summary(cdc$weight)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    68.0   140.0   165.0   169.7   190.0   500.0

summarise(
  cdc, 
  media_pe = mean(weight), 
  dp_pe = sd(weight), 
  media_alt = mean(height), 
  dp_alt = sd(height)
  ) %>%
  kable() %>%
  kable_styling(
    full_width = F, 
    bootstrap_options = c("striped", "hover", "condensed",
                          "responsive")
    )

mean(cdc$weight)

## [1] 169.683

median(cdc$height)

## [1] 67

cdc %>%
  group_by(gender) %>%
  summarise(
    media_pe = mean(weight), 
    dp_pe = sd(weight), 
    media_al = mean(height), 
    dp_al = sd(height)
    ) %>%
  kable() %>%
  kable_styling(
    full_width = F, 
    bootstrap_options = c("striped", "hover", "condensed",
                          "responsive")
    )

cdc %>%
  group_by(exerany) %>%
  summarise(
    median_pe = median(weight), 
    dq_peso = IQR(weight)
    ) %>%
  kable() %>%
     kable_styling(
    full_width = F, 
    bootstrap_options = c("striped", "hover", "condensed",
                          "responsive")
    )

cdc %>%
  select(smoke100) %>%
  table()

## smoke100
##     0     1 
## 10559  9441

cdc %>%
  select(smoke100) %>%
  table()/20000

## smoke100
##       0       1 
## 0.52795 0.47205

cdc %>%
  select(smoke100, gender) %>%
  mutate(smoke100 = ifelse(smoke100 == 1, "Sim", "Não")) %>%
  ggplot(aes(x = smoke100, fill = gender)) + 
  geom_bar(position = "dodge", width = 0.5) + 
  labs(
    x = "Fumou pelo menos 100 cigarros?",
    y = "Número de indivíduos",
    title = "Distribuição dos Fumantes por Gênero"
  ) + 
  theme(text = element_text(size = 18))

cdc %>%
  select(gender) %>%
  table() %>%
  prop.table()

## gender
##       m       f 
## 0.47845 0.52155

cdc %>%
  select(genhlth) %>%
  table() %>%
  prop.table()

## genhlth
## excellent very good      good      fair      poor 
##   0.23285   0.34860   0.28375   0.10095   0.03385

cdc %>%
  select(gender, smoke100) %>%
  table()

##       smoke100
## gender    0    1
##      m 4547 5022
##      f 6012 4419

cdc %>%
  select(gender, smoke100) %>%
  table() %>%
  prop.table(margin = 2)

##       smoke100
## gender         0         1
##      m 0.4306279 0.5319352
##      f 0.5693721 0.4680648

cdc %>%
  select(gender, smoke100) %>%
  table() %>%
  prop.table()

##       smoke100
## gender       0       1
##      m 0.22735 0.25110
##      f 0.30060 0.22095

nrow( cdc[cdc$smoke100 == 1 & cdc$gender == 'f', ] ) / nrow( cdc[cdc$smoke100 == 1, ] )

## [1] 0.4680648

cdc %>%
  select(gender, smoke100) %>%
  table() %>%
  prop.table(margin = 1)

##       smoke100
## gender         0         1
##      m 0.4751803 0.5248197
##      f 0.5763589 0.4236411

cdc %>%
  select(gender, smoke100) %>%
  table() %>%
  margin.table(margin = 1)

## gender
##     m     f 
##  9569 10431

cdc %>%
  select(gender, smoke100) %>%
  table() %>%
  margin.table(margin = 2)

## smoke100
##     0     1 
## 10559  9441

Visualizando Relações

cdc %>%
  ggplot(aes(x = height, y = weight)) + 
  geom_point(size = 1.5, color = "red") + 
  labs(
    title = "Peso versus altura"
  ) + 
  theme(text = element_text(size = 18))

cor(cdc$weight, cdc$height)

## [1] 0.5553222

cdc %>%
  select(smoke100, gender) %>%
  mutate(smoke100 = ifelse(smoke100 == 1, "Sim", "Não")) %>%
  ggplot(aes(x = smoke100, fill = gender)) + 
  geom_bar(position = "dodge", width = 0.5) + 
  coord_flip() + 
  labs(
    x = "Fumou pelo menos 100 cigarros?",
    y = "Número de indivíduos",
    title = "Distribuição dos Fumantes por Gênero"
  ) + 
  theme(text = element_text(size = 18))

cdc %>%
  ggplot(aes(x = gender, y = height, fill = gender)) + 
  geom_boxplot(width = 0.4) + 
  labs(
    x = "",
    title = "Altura em função do gênero"
  ) +
  theme_void() + 
  theme(text = element_text(size = 20))

cdc = mutate(cdc, imc = (weight/height^2)*703)

cdc %>%
  ggplot(aes(x = genhlth, y = imc, fill = genhlth)) + 
  geom_boxplot(size = 0.6) + 
  theme_void() + 
  scale_fill_manual(values = c("green", "skyblue", "orange", "yellow", "red")) + 
  theme(text = element_text(size = 20))

cdc %>%
  ggplot(aes(x = gender, y = imc, fill = gender)) + 
  geom_boxplot(size = 0.6, width = 0.4) + 
  theme_void() + 
  theme(text = element_text(size = 20))

Subconjuntos

É bastante comum agruparmos em subconjuntos elementos individuais do banco de dados que possuam características específicas. Nós podemos fazer isso facilmente por meio da função filter e uma série de operadores lógicos. Os operadores lógicos mais comuns na análise de dados são

• == - significa “igual a”
• != - significa “diferente de”
• > or < - significa “maior que” ou “menor que”
• >= or <= - siginifica “maior ou igual que” ou “menor ou igual que”

Usando esses operadores, nós podemos criar um subconjunto do banco de dados cdc relativo apenas aos homens e, então, salvar esta filtragem em um novo conjunto de dados chamado homens:

men = filter(cdc, gender == "m")

Aqui nós criamos um novo objeto, chamado men. Nós podemos dar uma olhada neste objeto clicando sobre o banco de dados homens que aparece em seu ambiente de traballho (Environment), ou podemos usar o comando head para dar uma espiada nas suas primeiras linhas:

head(men) %>%
  kable() %>%
  kable_styling(
    full_width = F, 
    bootstrap_options = c("striped", "hover", "condensed",
                          "responsive")
    )

Este novo conjunto de dados contém as mesmas variáveis de antes, apenas com um número menor de linhas. Nós também poderíamos dizer ao R para manter algumas variáveis específicas usando a função select que foi introduzida anteriormente, mas isto não é relevante no momento, uma vez que iremos focar na criação de subconjuntos do banco de dados.

men2 = cdc %>%
  subset(gender == "m", select = c(gender, genhlth, age))

Você poderá utilizar diversas dessas condições conjuntamente aos condicionadores & and |. O condicionador & corresponde à intersecção. A grosso modo, podemos dizer que & é equivalente a e. Nos comandos abaixo, agrupamos em um banco de dados os indivíduos que são homens e possuem mais de trinta anos:

homens_e_30mais = filter(cdc, gender == "m" & age > 30)

Já o condicionador | corresponde à união e, podemos dizer, que é equivalente a ou. O comando abaixo agrupa os indivíduos do conjunto de dados que são homens ou possuam mais de 30 anos:

homens_ou_30mais = filter(cdc, gender == "m" | age > 30)

À princípio, você pode usar esses condicionadores na quantidade e da maneira que decidir para formar subconjuntos. Se estes agrupamentos tem interpretações práticas é uma outra questão

• Crie um novo objeto menos23_fumante relativo àqueles que possuem menos de 23 anos e já fumaram 100 cigarros em seu tempo de vida. Escreva o comando que você utilizou para criar este novo objeto e apresente o número de casos que atendem a esse critério.

menos23_fumante = filter(cdc, age < 23 & smoke100 == 1)

Vemos então que 620 indivíduos do nosso banco de dados atendem à condição de ter menos que 23 anos e já ter fumado pelos menos 100 cigarros ao longo da vida.

Recapitulando

A esta altura, fizemos uma boa primeira exposição sobre a análise das informações contidas no questionário BRFSS. Nós descobrimos uma associação interessante entre fumo e gênero, e nós podemos comentar algo a respeito da relação entre a avaliação de saúde em geral dada pelas pessoas e seu próprio IMC. Nós também nos apropriamos de ferramentas computacionais essenciais – estatísticas descritivas, subconjuntos e gráficos – que nos servirão bem ao longo deste curso.

Por sua conta

1. Crie um gráfico de dispersão da variável peso em relação ao peso desejado. Descreva a relação entre essas duas variáveis.

cdc %>%
  ggplot(aes(x = wtdesire, y = weight)) + 
  geom_point(size = 1.8, color = "orange") +
  labs(
    title = "Peso desejado versus peso real",
    x = "peso real",
    y = "peso desejado"
  ) +
  theme(text = element_text(size = 18))

cor(cdc$wtdesire, cdc$weight)

## [1] 0.8000521

O gráfico de dispersão mostrado acima mostra uma relação praticamente linear entre as variáveis weight e wtdesire. Este fato é corroborado pelo coeficiente de correlação linear entre as variáveis, que foi de aproximadamente 0.8.

2. Vamos considerar uma nova variável: a diferença entre o peso desejado (wtdesire) e o peso atual (weight). Crie esta nova variável subtraindo as duas colunas na base de dados e atribuindo-as a um novo objeto chamado pdif.

pdif = cdc$wtdesire - cdc$weight

3. Que tipo de dados estão contidos na variável pdif? Se uma observação de pdif é 0, o que isso implica em relação ao peso atual e desejado de uma pessoa? E se o valor de pdif for positivo ou negativo?

typeof(pdif)

## [1] "integer"

typeof(cdc$weight)

## [1] "integer"

typeof(cdc$wtdesire)

## [1] "integer"

Como podemos ver na saída dos comandos acima, as variáveis weight e wtdesire são ambas do tipo integer, de modo que pdif, que é a diferença entre aquelas duas, é também do tipo inteiro. Vemos então que tanto o peso como o peso desejado foram registrados em quilos no nosso banco de dados, sem a parte fracionária.

Vemos ainda que, se pdif é zero, a pessoa já acredita que está em seu peso ideal, desejado. Se essa quantidade é positiva, ela acredita que seu peso precisa aumentar e, se é negativa, ela acredita que precisa perder peso.

4. Descreva a distribuição de pdif em termos de sua posição central, forma e variação, incluindo qualquer gráfico que você usar. O que isso nos diz sobre como as pessoas se sentem a respeito do seu peso atual?

cdc = mutate(cdc, pdif = wtdesire - weight)

cdc %>%
  select(pdif) %>%
  summary()

##       pdif        
##  Min.   :-300.00  
##  1st Qu.: -21.00  
##  Median : -10.00  
##  Mean   : -14.59  
##  3rd Qu.:   0.00  
##  Max.   : 500.00

cdc %>%
  ggplot(aes(x = pdif)) +
  geom_histogram(aes(y = after_stat(density)), 
                 fill = "orange", bins = 50, alpha = 0.8, color = "black") +
  labs(
    x = "P_dif", 
    y = "Densidade", 
    title = "Histograma da diferença entre o peso desejado e o atual"
  ) + 
  theme(text = element_text(size = 18))

Com este resumo dos dados da variável pdif que fizemos logo acima, percebemos que, em média, os participantes acreditam que estão próximos de seu peso ideal, com cerca de 6 quilos acima. Além disso, a partir do histograma e do boxplot gerados, podemos notar também que a distribuição das diferenças entre os pesos ideal e real dos participantes do estudo é praticamente simétrica. Contudo, podemos notar também que existem mais pessoas acreditando que estão acima do peso ideal do que abaixo dele.

5. Utilizando estatísticas numéricas e uma sequência de boxplots, determine se homens tendem a ver seu peso diferentemente das mulheres.

cdc %>%
  group_by(gender) %>%
  summarize(
    media_pdif = mean(pdif), 
    sd_pdif = sd(pdif), 
    mediana_pdif = median(pdif), 
    dq_pdif = IQR(pdif)
    ) %>%
    kable() %>%
  kable_styling(
    full_width = F, 
    bootstrap_options = c("striped", "hover", "condensed",
                          "responsive")
    )

cdc %>%
  ggplot(aes(x = gender, y = pdif, fill = gender)) +
  geom_boxplot(width = 0.4) + 
  theme_void() + 
  theme(text = element_text(size = 20))

Conforme vemos nas estatísticas calculadas acima, bem como nos boxplots, dentre os participantes do estudo, as mulheres tendem a achar que estão mais acima do peso do que os homens pensam sobre eles mesmos.

6. Agora chegou a hora de usar a criatividade. Encontre a média e o desvio padrão de peso e determine qual a proporção de pesos que estão a um desvio padrão da média. E a dois desvios padrões da média?

Para responder a essas perguntas, basta calcularmos, dentro do nosso dataframe, a proporção de elementos que satisfazem as seguintes condições

\[ |P-\bar{P}| < sd(P) \]

\[ |P-\bar{P}| < 2*sd(P) \]

nrow(filter(cdc, abs(weight - mean(weight)) < sd(weight)))/nrow(cdc)

## [1] 0.7076

nrow(filter(cdc, abs(weight - mean(weight)) < 2*sd(weight)))/nrow(cdc)

## [1] 0.9563

Como podemos ver, aproximadamente 70% dos valores dos pesos dos entrevistados estão dentro de 1 desvio padrão da média, ao passo que aproxidamente 95% dos valores estão dentro de dois desvios padrões de distância da média. Vemos, portanto, que o peso dos indivíduos entrevistados seguem uma distribuição aproximadamente normal. Isto pode ser corroborado com um histograma para essa variável:

cdc %>%
  ggplot(aes(x = weight)) +
  geom_histogram(
    aes(y = after_stat(density)), 
    fill = "blue",
    alpha = 0.5,
    col = "black"
    ) + 
  labs(
    x = "Peso", 
    y = "Densidade",
    title = "Histograma para os pesos dos participantes do estudo"
    ) + 
  theme(text = element_text(size = 18))

Conteúdos abordados neste laboratório

• Probabilidade
• Eventos independentes
• Simulações de eventos independentes no R

Mãos Quentes (Hot Hands)

Jogadores de basquete que pontuam várias vezes seguidas costumam ser descritos como tendo as mãos quentes. Fãs e jogadores acreditam há muito tempo no fenômeno das mãos quentes, que refuta o pressuposto de que cada lance é independente do próximo. Contudo, um artigo de 1985 escrito por Gilovich, Vallone e Tversky coletou evidências que contradizem essa crença e mostrou que lances sucessivos são eventos independentes. Este artigo iniciou uma grande controvérsia que continua até hoje, como você pode verificar se procurar por hot hand basketball no Google.

Não temos a expectativa de resolver esta controvérsia hoje. Entretanto, neste laboratório nós aplicaremos um procedimento para responder a questões como essa. Os objetivos deste laboratório são (1) refletir sobre o efeito de eventos dependentes e independentes, (2) aprender como simular sequências de lances no R, e (3) comparar a simulação com os dados efetivos para determinar se o fenômeno das mãos quentes parece ser real.

Começando

Nossa investigação terá como foco a performance de um jogador: Kobe Bryant do Los Angeles Lakers. Sua performance contra o Orlando Magic nas finais de 2009 da NBA lhe deram o título de Jogador Mais Valioso (ou, em inglês, Most Valuable Player) e vários espectadores comentaram como ele parecia demonstrar a característica de mãos quentes. Vamos carregar alguns dados desses jogos que iremos utilizar durante este laboratório.

load(file = "hot_hand.Rdata")

Você deve ver dois objetos na aba Environment: kobe_basket e calc_streak:

• kobe_basket: um conjunto de dados de 133 observações e 6 variávies, onde cada linha registra um arremesso feito por Kobe Bryant. A variável basket indica se o arremesso foi um acerto (H) ou um erro (M).

• calc_streak: uma função customizada que irá nos ajudar a calcular comprimentos de sequências de acertos (mais sobre isso à frente).

names(kobe_basket)

## [1] "vs"          "game"        "quarter"     "time"        "description"
## [6] "shot"

head(kobe_basket) %>%
  kable() %>%
  kable_styling(
    full_width = F, 
    bootstrap_options = c("striped", "hover", "condensed",
                          "responsive")
  )

Apenas olhando para a sequência de acertos e erros pode ser difícil de auferir se é provável que Kobe estava arremessando com as mãos quentes. Uma possível maneira de abordarmos este problema é supormos que arremessadores com as mãos quentes tendem a conseguir uma longa sequência de acertos. Para este laboratório, definiremos o comprimento de uma sequência de acertos como o número de cestas consecutivas até acontecer um erro.

Por exemplo, no Jogo 1, Kobe teve a seguinte sequência de acertos e erros de suas nove tentativas de arremessos no primeiro quarto:

Para verificar esta informação, use o seguinte comando:

kobe_basket%>%
  filter(game == 1, quarter == 1) %>%
  select(shot)

##   shot
## 1    H
## 2    M
## 3    M
## 4    H
## 5    H
## 6    M
## 7    M
## 8    M
## 9    M

Dentre as nove tentativas de arremesso existem seis sequências, que são separadas por um “|” acima. Seus comprimentos são um, zero, dois, zero, zero, zero (em ordem de ocorrência). Vemos então que em uma sequência de comprimento 1 existe 1 acerto e 1 erro, ao passo que na sequência de comprimento 0 existe 0 acerto e 1 erro. Uma sequência de comprimento 2, por sua vez, contêm 2 acertos e 1 erro. Portanto, uma sequência de comprimento L contêm L acertos e 1 erro.

Contar os comprimentos de todos os 133 arremessos manualmente seria extremamente tedioso, então iremos utilizar a função calc_streak para realizar este trabalho para nós. Os resultados serão armazenados em um conjunto de dados chamado kobe_streak como a variável length (comprimento em inglês).

kobe_streak = kobe_basket %>%
  select(shot) %>%
  calc_streak()

Nós podemos dar uma olhada na distribuição desses comprimentos com o comando abaixo.

kobe_streak %>%
  ggplot(aes(x = length)) +
  geom_histogram(binwidth = 1, col = "black", fill = "orange") + 
  labs(
    title = "Tamanho das sequências de arremesso de Kobe",
    y = "N"
    ) +
  theme(text = element_text(size=18))

summary(kobe_streak)

##      length      
##  Min.   :0.0000  
##  1st Qu.:0.0000  
##  Median :0.0000  
##  Mean   :0.7632  
##  3rd Qu.:1.0000  
##  Max.   :4.0000

table(kobe_streak)

## length
##  0  1  2  3  4 
## 39 24  6  6  1

Conforme podemos ver na saída de código acima, o comprimento típico das sequências de Kobe foi zero. Além disso, podemos ver também que a mediana dos comprimentos das sequências é também zero, enquanto a média é um número menor que 1.

Ao que comparar

Nós mostramos que Kobe teve algumas sequências longas de arremesso, mas elas são longas o suficiente para apoiar a crença de que ele tinha mãos quentes? Com o que podemos compará-las?

Para responder a essa pergunta, vamos retornar ao conceito de independência. Dois eventos são independentes se o resultado de um evento não afeta o resultado do outro. Se cada arremesso do jogador for um evento independente, acertar ou errar o primeiro arremesso não afetará a probabilidade de ele converter ou errar o segundo arremesso.

Um arremessador com as mãos quentes terá arremessos que não são independentes um do outro. Mais especificamente, se o arremessador converte seu primeiro arremesso, o modelo acerca das mãos quentes afirma que ele terá uma probabilidade maior de converter seu segundo arremesso.

Vamos supor por um momento que o modelo das mãos quentes é válido para Kobe. Durante sua carreira, o percentual de vezes que Kobe acertou uma cesta (ou seja, seu percentual de arremessos convertidos) é de cerca de 45%, ou, em notação de probabilidade,

Se ele converte o primeiro arremesso e tem as mãos quentes (arremessos não independentes), então a probabilidade de ele converter seu segundo arremesso deveria aumentar para, digamos, 60%,

Como um resultado do aumento da probabilidade, seria esperado que Kobe obtivesse sequências mais longas. Compare com a perspectiva cética de que Kobe não tenha as mãos quentes, ou seja, que cada arremesso é independente do anterior. Se ele acerta seu primeiro arremesso, a probabilidade de ele acertar o segundo continua sendo 0.45.

Em outras palavras, converter o primeiro arremesso não afeta de maneira alguma a probabilidade de ele converter o segundo arremesso. Se os arremessos de Kobe forem independentes, então ele teria a mesma probabilidade de acertar cada arremesso independentemente de seus arremessos anteriores: 45%.

Agora que reformulamos a situação em termos de arremessos independentes, vamos retornar à questão: como podemos saber se as sequências de arremessos de Kobe são longas o suficiente para indicar que ele tem mãos quentes? Podemos comparar o tamanho de suas sequências a alguém que não tem as mãos quentes: um arremessador independente.

Simulações no R

Embora não tenhamos dados de um “arremessador independente”, esse tipo de dado é muito fácil de ser simulado no R. Numa simulação, você define as regras básicas de um processo aleatório e então o computador utiliza números aleatórios para gerar um resultado fiel a essas regras. Como um exemplo simples, vamos simular o lançamento de uma moeda honesta onde K simboliza cara e C simboliza coroa:

face_moeda = c("K", "C")
sample(face_moeda, size = 1, replace = T)

## [1] "K"

Primeiramente registramos os possíveis resultados no vetor face_moeda. O comando sample realiza uma amostra do experimento ao selecionar aleatoriamente um dos valores K ou C. O vetor face_moeda pode ser entendido como um chapéu com duas tiras de papel dentro dele: numa tira está escrito K e na outra C. A função sample sorteia uma tira de dentro do chapéu e revela se ela é cara ou coroa.

Execute o segundo comando listado acima várias vezes. Da mesma maneira que quando jogamos uma moeda, algumas vezes você obterá cara, algumas vezes você obterá coroa, mas a longo prazo você esperaria obter um número mais ou menos igual de cada possível resultado.

Se você quisesse simular o lançamento de uma moeda honesta 100 vezes, você poderia rodar a função 100 vezes ou, de forma mais simples, ajustar o argumento size (tamanho), que regula quantas amostras retirar. O argumento replace indica se há (TRUE) ou não (FALSE) reposição dos possíveis valores. Voltando à analogia do chapéu, TRUE indica que nós recolocamos a tira de papel de volta no chapéu antes de retirar as amostras seguintes. Salve o vetor resultante de caras e coroas num novo objeto denominado sim_moeda_honesta.

sim_moeda_honesta = sample(face_moeda, 100, replace = T)

sim_moeda_honesta

##   [1] "C" "C" "K" "K" "K" "K" "C" "K" "K" "K" "K" "K" "K" "K" "C" "K" "K" "C"
##  [19] "C" "C" "K" "K" "C" "K" "C" "K" "K" "K" "K" "K" "K" "K" "K" "C" "C" "C"
##  [37] "K" "K" "C" "C" "C" "C" "C" "C" "K" "K" "K" "C" "K" "K" "K" "C" "K" "K"
##  [55] "C" "C" "K" "C" "K" "K" "K" "C" "C" "K" "C" "C" "C" "C" "K" "K" "K" "K"
##  [73] "K" "C" "C" "C" "K" "K" "K" "C" "K" "C" "K" "K" "C" "K" "K" "C" "K" "C"
##  [91] "K" "C" "C" "K" "C" "K" "K" "K" "C" "C"

Para visualizar os resultados desta simulação, digite o nome do objeto e então use o comando table pra contar o número de caras e coroas.

table(sim_moeda_honesta)

## sim_moeda_honesta
##  C  K 
## 42 58

Uma vez que há apenas dois elementos no vetor face_moeda, a probabilidade de que um lançamento dê o resultado cara é 0,5. Digamos que estamos tentando simular uma moeda viciada, onde a probabilidade de sair cara é 20%. Podemos ajustar a amostra adicionando um argumento denominado prob, que fornece um vetor com as respectivas probabilidades.

sim_moeda_viciada = sample(face_moeda, size = 100, replace = T, prob = c(0.2, 0.8))

prob = c(0.2, 0.8) indica que, para os dois elementos no vetor face_moeda, nós queremos selecionar o primeiro (cara) com probabilidade 0,2, e o segundo (coroa), com probabilidade 0,8. Outra maneira de pensar sobre esse cenário é imaginar o espaço amostral como um saco contendo 10 fichas, sendo 2 marcadas como “cara” e 8 como “coroa”. Portanto, a cada seleção, a probabilidade de retirar uma ficha escrito “cara” é 20%, e “coroa” é 80%.

Uma nota sobre sementes: Especificar uma semente irá fazer com que o R obtenha a mesma amostra cada vez que você aperta o botão Knit. Isso garante a reprodutibilidade do seu trabalho, uma vez que os seus resultados serão os mesmos cada vez que você compilar os seus comandos. Você pode especificar uma semente com o comando abaixo:

set.seed(123) #Podemos usar qualquer valor como semente, desde que 
#utilizemos sempre o mesmo ao longo das nossas simulações.

set.seed(123)

sim_moeda_viciada_reprod = sample(face_moeda, 100, replace = T, prob = c(0.2,0.8))

sim_moeda_viciada_reprod

##   [1] "C" "C" "C" "K" "K" "C" "C" "K" "C" "C" "K" "C" "C" "C" "C" "K" "C" "C"
##  [19] "C" "K" "K" "C" "C" "K" "C" "C" "C" "C" "C" "C" "K" "K" "C" "C" "C" "C"
##  [37] "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "K" "C" "C" "C" "C"
##  [55] "C" "C" "C" "C" "K" "C" "C" "C" "C" "C" "K" "C" "K" "K" "C" "C" "C" "C"
##  [73] "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "K" "K" "K" "C"
##  [91] "C" "C" "C" "C" "C" "C" "C" "C" "C" "C"

table(sim_moeda_viciada_reprod)

## sim_moeda_viciada_reprod
##  C  K 
## 82 18

Num certo sentido, nós reduzimos o tamanho da tira de papel que diz “cara”, tornando-o menos provável de ser escolhido, e nós aumentamos o tamanho da tira de papel que diz “coroa”, tornando-o mais provável de ser retirado. Quando simulamos a moeda honesta, ambas as tiras de papel tinham o mesmo tamanho. Isso acontece de maneira padrão se você não fornecer o argumento prob; todos os elementos no vetor face_moeda terão igual probabilidade de serem escolhidos.

Simulando o arremessador independente

Para simular um jogador de basquete que arremessa de forma independente, utilizamos o mesmo mecanismo que empregamos para simular o lançamento de uma moeda. Para simular um único arremesso de um arremessador independente, com um percentual de acerto de 50%, digitamos

arremesso = c("H", "M")
sim_cesta = sample(arremesso, 1, replace = T)

Para fazermos uma comparação válida entre Kobe e nosso arremessador independente simulado, precisamos alinhar tanto seus percentuais de arremesso quanto seus números de arremessos tentados.

Considerando os 133 arremessos de Kobe nas finais de 2009 da NBA, bem como seu percentual histório de acertos, que era de 45%, podemos simular como seria sua versão independente da seguinte forma:

set.seed(123)
sim_cesta = sample(arremesso, size = 133, replace = T, prob = c(0.45, 0.55))

Observe que nomeamos o novo vetor como sim_cesta, com o mesmo nome que demos ao vetor anterior correspondente a um percentual de arremesso de 50%. Nessa situação, o R sobrescreve o objeto antigo com o novo, portanto sempre se certifique que você não precisa da informação no vetor antigo antes de atribuir um novo objeto ao seu nome.

Com os resultados da simulação salvos como sim_cesta, temos os dados necessários para comparar Kobe ao nosso arremessador independente. Podemos visualizar os dados de Kobe em conjunto com os nossos dados simulados.

Ambos os conjuntos de dados representam o resultado de 133 tentativas de arremessos, cada uma com o mesmo percentual de acerto de 45%. Sabemos que nossos dados simulados são de um arremessador que arremessa de forma independente. Quer dizer, sabemos que o arremessador simulado não tem as mãos quentes.

Por sua conta

1. Utilizando a função calc_streak, calcule os comprimentos das sequências de arremessos convertidos de sim_cesta e salve os resultados em um conjunto de dados chamado sim_streak. Note que sim_streak é apenas um vetor e não uma variável de um conjunto de dados como o que utilizamos para calcular as sequências de Kobe. Desse modo, não é preciso selecionar a variável dentro de um conjunto de dados, bastando apenas aplicar a função ao vetor.

sim_streak = calc_streak(sim_cesta)

2. Descreva a distribuição das sequências de arremessos. Qual é o comprimento de sequência típico para o arremessador independente simulado com um percentual de arremesso de 45%? Qual o comprimento da sequência mais longa de cestas em 133 arremessos? Não esqueça de incluir um gráfico à sua resposta.

Para começar, vamos ver um resumo estatístico das sequências de Kobe e do arremessador independente:

summary(kobe_streak)

##      length      
##  Min.   :0.0000  
##  1st Qu.:0.0000  
##  Median :0.0000  
##  Mean   :0.7632  
##  3rd Qu.:1.0000  
##  Max.   :4.0000

summary(sim_streak)

##      length      
##  Min.   :0.0000  
##  1st Qu.:0.0000  
##  Median :0.0000  
##  Mean   :0.8108  
##  3rd Qu.:1.0000  
##  Max.   :7.0000

Vemos então que o tamanho máximo de uma sequência de Kobe foi igual a 4, ao passo que o arremessador independente alcançou um máximo de 7. As medianas são iguais e as médias muito parecidas.

Podemos ainda executar o comando table, combinado com o comando prop.table(), a fim de verificar, em cada conjunto de dados, o tamanho de sequência que mais ocorre:

prop.table(table(kobe_streak))

## length
##          0          1          2          3          4 
## 0.51315789 0.31578947 0.07894737 0.07894737 0.01315789

prop.table(table(sim_streak))

## length
##          0          1          2          3          4          7 
## 0.54054054 0.28378378 0.08108108 0.05405405 0.02702703 0.01351351

Vemos, portanto, que as sequências ocorrem com frequências muito parecidas nos dois conjuntos de dados. Por exemplo, aproximadamente 51% das sequências de Kobe naquela final de temporada teve comprimento zero. No caso do simulador independente, esse percentual foi perto de 54%. Já as sequências de tamanho 2 ocorreram com frequência de aproximadamente 8% em ambos os casos. Por meio dessa análise, podemos concluir que os dois arremessadores são compatíveis.

Por fim, vamos comparar os histogramas para os comprimentos dos arremessos de Kobe e do arremessador independente:

hist_kobe = ggplot(kobe_streak, aes(x = length)) +
  geom_histogram(col = "black", fill = "green", binwidth = 1) +
  labs(
    x = "Comprimento",
    y = "N",
    title = "Sequências de Kobe"
  ) +
  theme(text = element_text(size=18))

hist_sim = ggplot(sim_streak, aes(x = length)) +
  geom_histogram(col = "black", fill = "orange", binwidth = 1) +
  labs(
    x = "Comprimento",
    y = "",
    title = "Sequências simuladas"
  ) +
  theme(text = element_text(size=18))

cowplot::plot_grid(hist_kobe, hist_sim, labels = "AUTO")

Conteúdos abordados neste laboratório

• Amostragem no R
• Distribuição amostral da média
• Distribuição normal
• Teorema do Limite Central
• Histogramas
• Box plots

Especificando uma semente: Nós iremos tomar amostras aleatórias e gerar algumas distribuições amostrais neste laboratório. Desse modo, você deverá especificar uma semente logo no início de seu relatório. Caso tenha dúvidas sobre como definir esta semente, retorne ao laboratório anterior onde explicamos este conceito.

set.seed(123)

Os dados

Vamos analisar dados do setor imobiliário da cidade de Ames, no estado de Iowa, Estados Unidos. Os detalhes de cada transação imobiliária na cidade de Ames são registrados pelo escritório da Secretaria Municipal da Receita da cidade. Nosso foco particular para este laboratório será todas as vendas de casa em Ames entre 2006 e 2010. Essa coleção representa nossa população de interesse. Neste laboratório queremos aprender sobre essas vendas de casas retirando pequenas amostras da população completa. Vamos importar os dados.

load("ames.rda")
names(ames)

##  [1] "Order"           "PID"             "area"            "price"          
##  [5] "MS.SubClass"     "MS.Zoning"       "Lot.Frontage"    "Lot.Area"       
##  [9] "Street"          "Alley"           "Lot.Shape"       "Land.Contour"   
## [13] "Utilities"       "Lot.Config"      "Land.Slope"      "Neighborhood"   
## [17] "Condition.1"     "Condition.2"     "Bldg.Type"       "House.Style"    
## [21] "Overall.Qual"    "Overall.Cond"    "Year.Built"      "Year.Remod.Add" 
## [25] "Roof.Style"      "Roof.Matl"       "Exterior.1st"    "Exterior.2nd"   
## [29] "Mas.Vnr.Type"    "Mas.Vnr.Area"    "Exter.Qual"      "Exter.Cond"     
## [33] "Foundation"      "Bsmt.Qual"       "Bsmt.Cond"       "Bsmt.Exposure"  
## [37] "BsmtFin.Type.1"  "BsmtFin.SF.1"    "BsmtFin.Type.2"  "BsmtFin.SF.2"   
## [41] "Bsmt.Unf.SF"     "Total.Bsmt.SF"   "Heating"         "Heating.QC"     
## [45] "Central.Air"     "Electrical"      "X1st.Flr.SF"     "X2nd.Flr.SF"    
## [49] "Low.Qual.Fin.SF" "Bsmt.Full.Bath"  "Bsmt.Half.Bath"  "Full.Bath"      
## [53] "Half.Bath"       "Bedroom.AbvGr"   "Kitchen.AbvGr"   "Kitchen.Qual"   
## [57] "TotRms.AbvGrd"   "Functional"      "Fireplaces"      "Fireplace.Qu"   
## [61] "Garage.Type"     "Garage.Yr.Blt"   "Garage.Finish"   "Garage.Cars"    
## [65] "Garage.Area"     "Garage.Qual"     "Garage.Cond"     "Paved.Drive"    
## [69] "Wood.Deck.SF"    "Open.Porch.SF"   "Enclosed.Porch"  "X3Ssn.Porch"    
## [73] "Screen.Porch"    "Pool.Area"       "Pool.QC"         "Fence"          
## [77] "Misc.Feature"    "Misc.Val"        "Mo.Sold"         "Yr.Sold"        
## [81] "Sale.Type"       "Sale.Condition"

Como você pode observar, o banco de dados está em sua língua originária, inglês, e existe uma quantidade consideravelmente grande de variáveis. Para este laboratório, restringiremos nossa atenção para somente duas variáveis: a área habitável da casa acima do nível do solo em pés quadrados (area) e o preço da venda (price) em dólares.

Nós podemos explorar a distribuição da área das casas em nossa população de vendas de casas calculando algumas estatísticas sumárias e criando um histograma. Vamos iniciar criando uma visualização por meio do histograma:

ames %>%
  ggplot(aes(x = area)) + 
  geom_histogram(fill = "orange", col = "black", binwidth = 250) + 
  labs(
    x = "Area",
    y = "N",
    title = "Histograma para as áreas das casas em Ames, Iowa"
  ) + 
  theme(text = element_text(size = 18))

Vamos obter, também, algumas medidas descritivas. Observe que podemos fazer isso por meio da função summarise. O número de estatísticas a serem calculadas fica a nosso critério. Algumas das funções abaixo já foram discutidas no primeiro laboratório (como mean, median, sd, IQR, min, e max). A função quantile aparece pela primeira vez em nossos laboratórios, correspondendo ao quantil especificado da distribuição de interesse. Por exemplo, a expressão quantile(x, 0.25) irá nos revelar o valor que acumula 25% (Q1) da distribuição da variável X. Conforme pode ser acompanhado em qualquer livro-texto de estatística básica, esses valores são úteis na descrição da distribuição dos dados: como exemplo, poderíamos estar interessados na “região central que acumula 50% das casas com área dada em pés quadrados”. Na saída abaixo, usamos também a função kable, do pacote kableExtra, para apresentar a saída da função summarise na forma de uma tabela:

n = nrow(ames)
ames %>%
  summarise(
    mu = mean(area), 
    pop_med = median(area),
    sigma = ((n-1)/n)*sd(area), 
    pop_iqr = IQR(area),
    pop_min = min(area), pop_max = max(area),
    pop_q1 = quantile(area, 0.25),# primeiro quartil, 25º percentil
    pop_q3 = quantile(area, 0.75) # terceiro quartil, 75º percentil
    ) %>%
  kable() %>%
  kable_styling(
    full_width = F, 
    bootstrap_options = c("striped", "hover", "condensed",
                          "responsive")
  )

• Exercício 1: Descreva a distribuição desta população usando visualizações e estatísticas descritivas. Você não precisa usar todas as estatísticas descritivas existentes, mas deve decidir quais delas são relevantes baseando-se na forma da distribuição. Não esqueça de incluir os gráficos e saídas das estatísticas calculadas no seu relatório juntamente com a sua narrativa.

ames %>%
  select(area) %>%
  filter(abs(area - mean(area)) < ((n-1)/n)*sd(area)) %>%
  nrow()/nrow(ames)

## [1] 0.7167235

Como podemos notar a partir do histograma plotado e da saída do código logo acima, assim como através das estatísticas descritivas selecionadas na função summarise, a distribuição das áreas das habitações do nosso conjunto de dados é aproximadamente normal, porém, com uma leve assimetria à direita, devido aos poucos imóveis que possuem área muito maior que os demais (outliers). Isso pode ser reforçado pela visualização de um boxplot:

ames %>%
  ggplot(aes(x = "", y = area)) + 
  geom_boxplot(width = 0.2, fill = "orange", outlier.color = "red") + 
  labs(
    title = "Boxplot para as áreas das casas em Ames, Iowa",
    x = "Area"
  ) + 
  theme_void() +
  theme(text = element_text(size = 18))

A distribuição Amostral Desconhecida

Neste laboratório nós temos acesso à população inteira, mas isso raramente acontece na vida real. Reunir informação sobre uma população inteira costuma ser muito custoso ou impossível. Por essa razão, frequentemente retiramos uma amostra da população e a utilizamos para compreender propriedades da população.

Se estivermos interessados em estimar a área habitável média na cidade de Ames com base numa amostra, podemos utilizar o seguinte comando para sondar a população.

set.seed(123)
amostra1 = slice_sample(ames, n = 50, replace = T)

Esse comando retira uma amostra aleatória simples de tamanho 50 do dataframe ames, que é atribuída à variável amostra1. É como se fossemos ao banco de dados da Secretaria Municipal da Fazenda e retirássemos os arquivos de 50 vendas de casas aleatoriamente. Trabalhar com esses 50 arquivos seria consideravelmente mais simples do que lidar com todas as 2930 vendas de casas.

• Exercício 2: Descreva a distribuição desta amostra. Como ela se compara à distribuição da população? Dica A função slice_sample retira uma amostra das observações (i.e. linhas) no banco de dados. Você pode utilizar os mesmos nomes para as variáveis do banco de dados. Os códigos utilizados no exercício anterior podem ser muito úteis para visualizar e descrever a distribuição, mas tenha cuidado ao nomear as estatísticas uma vez que mu e sigma são parâmetros populacionais já definidos e o que você tem é uma amostra e, consequentemente, estatísticas amostrais. Você pode nomear à sua escolha as estatísticas calculadas, sendo muito útil indicar no nome atribuído que seus valores são oriundos de uma amostra.

Se estamos interessados em estimar a área habitável média nas casas da cidade de Ames utilizando esta amostra, nosso melhor palpite é a média da amostra.

summarise(amostra1, x_bar = mean(area))

## # A tibble: 1 × 1
##   x_bar
##   <dbl>
## 1 1459.

Dependendo de quais foram as 50 casas que foram sorteadas, sua estimativa pode estar um pouco acima ou abaixo da média populacional verdadeira de 1499,69 pés quadrados. De maneira geral, mesmo assim, a média da amostra costuma ser uma estimativa muito boa da média da área habitável, e nós a obtemos por meio de uma amostra de menos de 3% da população.

Não é surpreendente que, a cada vez que retiramos uma nova amostra aleatória, obtemos uma média amostral diferente. É útil ter uma ideia de quanta variabilidade podemos esperar quando estimamos a média populacional desta maneira. A distribuição das médias amostrais, denominada de distribuição amostral, pode nos ajudar a compreeder essa variabilidade. Neste laboratório, uma vez que temos acesso à população, podemos elaborar a distribuição amostral para a média amostral repetindo os passos acima várias vezes. Agora iremos gerar 2000 amostras e calcularemos a média amostral de cada uma delas. Note que utilizamos o argumento replace = TRUE pelo simples fato de que amostras aleatórias simples são obtidas da distribuição original para cada réplica.

set.seed(123)
amostra50 = rdply(2000, slice_sample(ames, n = 50, replace = T), .id = "id_sample")

para determinar as médias de cada uma das 2000 amostras de 50 elementos, utilizamos as funções da biblioteca dplyr, além do operador %>%:

medias_amostra50 = amostra50 %>%
  group_by(id_sample) %>%
  summarise(x_bar = mean(area))

Interlúdio: A função rdply

Na seção anterior, utilizamos a função rdply, do pacote plyr, para gerar L = 2000 amostras de tamanho n = 50 do nosso conjunto de dados ames. Para tanto, continuamos a utilizar a função slice_sample, da biblioteca dplyr. Contudo, como nesse caso não precisávemos apenas de uma amostra, mas de 2000, selecionadas aleatoriamente do nosso dataframe, foi necessário a utilização da função rdply.

A função rdplypermite replicar o cálculo de uma expressão sem a necessidade de utilizar uma estrutura de repetição (ou loop), de modo a tornar a programação mais limpa e rápida. Além disso, os resultados dos cálculos efetuados são alocados em um novo dataframe, associados a um identificador (argumento .id, inserido na função rdply).

Para separar o dataframe criado em 2000 amostras isoladas, podemos usar a função split, da biblioteca plyr. Poderíamos fazer, por exemplo, amostra50 = split(amostra50, rep(1:2000, each = nrow(amostra50)/2000)). Neste caso, seria criada uma lista com os 2000 dataframes. Essa talvez não seja a opção mais interessante, uma vez que irá ocupar muito espaço de memória. Uma alternativa, por exemplo, é fazer o slice do dataframe que é retornado da função rdply, por meio da coluna .id.

Tamanho da Amostra e a Distribuição Amostral

Deixando de lado os aspectos técnicos de programação, vamos retomar a razão pela qual utilizamos a função rdply: calcular uma distribuição amostral, especificamente, esta aqui.

medias_amostra50 %>%
  ggplot(aes(x = x_bar)) + 
  geom_histogram(
    aes(y = after_stat(density)), 
    binwidth = 10, 
    fill = "green", 
    col = "black") + 
  stat_function(
    fun = dnorm,
    args = 
      list(mean = mean(medias_amostra50$x_bar), sd = sd(medias_amostra50$x_bar)), 
    linewidth = 1, 
    color = "red"
    ) + 
  theme(text = element_text(size = 18)) + 
  labs(
    title = "Distribuição das médias amostrais",
    subtitle = "2000 amostras com n = 50 cada"
  )

A distribuição amostral que calculamos nos informa bastante sobre as estimativas da área habitável das casas na cidade de Ames. Uma vez que a média amostral é um estimador não viesado, a distribuição amostral estará centrada na verdadeira média da área habitável da população, e a dispersão da distribuição indica quanta variabilidade é possível ao se amostrar somente 50 vendas de casas.

Para ter uma ideia melhor do efeito do tamanho da amostra na distribuição amostral, vamos construir distribuições a partir de amostras de diferentes tamanhos:

Passo 1: Crie três distribuições amostrais, cada uma contendo 2000 médias amostrais, oriundas de amostras de tamanho \(n=10\), \(n=50\) e \(n=100\).

set.seed(123)
medias_amostra10 = rdply(2000, slice_sample(ames, n = 10, replace = T), .id = "id_sample") %>%
  group_by(id_sample) %>%
  summarise(x_bar = mean(area))

set.seed(123)
medias_amostra100 = rdply(2000, slice_sample(ames, n = 100, replace = T), .id = "id_sample") %>%
  group_by(id_sample) %>%
  summarise(x_bar = mean(area))

Passo 2: Construa o gráfico dessas três distribuições amostrais, uma seguida da outra e, na mesma escala, de modo que possamos comparar as formas, centros e dispersões de cada uma.

• Passo 2a: Combine essas três distribuições amostrais (cada uma com 2000 valores) em um novo conjunto de dados contendo 6000 valores. Para combinar bancos de dados, linha por linha, nós usamos a função rbind.

• Passo 2b: Adicione uma nova coluna chamada tamanho_amostra aos dados que você acabou de criar, indicando o tamanho da amostra utilizado ao obter cada entrada. Lembre-se que os primeiros 2000 valores foram obtidos de amostras de tamanho 10, os 2000 seguintes de amostras de tamanho 50 e, finalmente, os 2000 últimos valores foram obtidos de amostras de tamanho 100. Portanto, essa nova variável é simplesmente a repetição do valor 10 por 2000 vezes, seguida de 50 por 2000 vezes e, novamente, 100 por 2000 vezes. O uso da função factor indica ao R que a variável é qualitativa, e não numérica. Lembre-se também de usar a função mutate para adicionar a nova variável ao banco de dados. O argumento rep(c(10,50,100), each=2000) diz ao R que os valores 10, 50 e 100 serão repetidos cada um 2000 vezes.

medias_amostra = rbind(medias_amostra10, medias_amostra50, medias_amostra100) %>%
  select(-id_sample) %>%
  mutate(tamanho_amostra = factor(rep(c(10,50,100), each = 2000)))

• Passo 2c: Finalmente, construa os três gráficos representando estas três distribuições amostrais. Podemos fazer isso via faceting, i.e. criar faces separadas no gráfico para cada uma das distribuições. Lembre-se que nós identificamos cada distribuição a partir da variável tamanho_amostra, que acabamos de criar. A expressão ~tamanho_amostra é nova: ela separa o eixo-y do nosso gráfico em três faces segundo o tamanho_amostra e não separa o eixo-x, o que significa que os gráficos serão dispostos um acima do outro ao invés de um ao lado do outro. Decidimos por esta abordagem para facilitar a comparação entre a posição e a dispersão das distribuições.

medias_amostra %>%
  ggplot(aes(x = x_bar)) + 
  geom_histogram(binwidth = 20, fill = "orange", col = "black") + 
  facet_wrap(~tamanho_amostra, ncol = 1, strip.position =  "right") +
  theme(text = element_text(size = 18))

Vemos, portanto, que o aumento do tamanho da amostra tem pouca influência sobre o centro da distribuição (a média amostral, que é uma estimativa da média populacional). A dispersão, porém, diminui bastante com o aumento do tamanho da amostra.

Isso mostra, na prática, o que diz o teorema do limite central: a distribuição amostral da média é normal com média igual à média populacional e variância igual a \(\sigma/\sqrt{n}\), em que \(n\) é o tamanho da amostra e \(\sigma\) é o desvio padrão populacional(MORETTIN; BUSSAB, 2017). Em outras palavras, podemos dizer que, se sorteamos aleatoriamente um número grande de amostras suficientemente grande, a distribuição da média amostral será normal, centrada na média da população e com variância inversamente proporcional ao tamanho das amostras.

Por Sua Conta

Até agora, nós nos ocupamos em estimar a média da área habitável nas casas do município de Ames. Agora você tentará estimar a média dos preços das casas.

1. Retire uma amostra aleatória de tamanho 15 da variável price e calcule a média. Com essa amostra, qual é sua melhor estimativa pontual para a média populacional?

Como sabemos, um bom estimador para a média populacional é a média amostral \(\bar{X}\), que é um estimador não viesado. Para esse caso, a média amostral seria:

set.seed(123)
ames %>%
  select(price) %>%
  slice_sample(n = 15, replace = T) %>%
  summarise(x_bar = mean(price))

## # A tibble: 1 × 1
##     x_bar
##     <dbl>
## 1 174223.

2. Já que você tem acesso à população, simule a distribuição amostral de \(\bar{X}_{price}\), retirando 2000 amostras de tamanho 15 da população e calculando 2000 médias amostrais. Armazene essas médias em um vetor com o nome medias_amostra15. Crie um gráfico com os resultados, e então descreva a forma dessa distribuição amostral. Baseado nessa distribuição amostral, qual seria seu palpite para a média dos preços das casas na população? Por fim, calcule e informe a média populacional.

set.seed(123)
medias_amostra15 = rdply(2000, slice_sample(ames, n = 15, replace = T), .id = "id_sample") %>%
  group_by(id_sample) %>%
  summarise(x_bar = mean(price))

medias_amostra15 %>%
  ggplot(aes(x = x_bar)) +
  geom_histogram(
    aes(y = after_stat(density)), 
    binwidth = 2000, 
    fill = "orange", 
    col = "black"
    ) + 
  labs(
    title = "Distribuição das médias amostrais dos preços",
    subtitle = "2000 amostras com n = 15 cada"
  ) + 
  theme(text = element_text(size = 18)) + 
  stat_function(
    fun = dnorm, 
    args = list(mean = mean(medias_amostra15$x_bar), sd = sd(medias_amostra15$x_bar)),
    linewidth = 1, 
    color = "red"
    )

Conforme podemos ver no histograma plotado logo acima, a forma da distribuição da média para essas 2000 amostras de tamanho 15 é de uma normal, como seria de esperar pelo teorema do limite central. Segundo esse histograma, a melhor estimativa para a média populacional seria em torno de $ 180.000,00.

Além disso, a média populacional dos preços dos imóveis do nosso conjunto de dados é igual a 180.796,06 dólares.

3. Mude o tamanho da sua amostra de 15 para 150, e então calcule a distribuição amostral utilizando o mesmo método descrito acima, e guarde as médias em um novo vetor com o nome medias_amostra150. Descreva a forma dessa distribuição amostral e compare-a com a distribuição amostral para a amostra de 15 elementos. Com base nessa distribuição amostral, qual seria seu palpite sobre a média dos preços de vendas de casas no município de Ames?

set.seed(123)
medias_amostra150 = rdply(2000, slice_sample(ames, n = 150, replace = T), .id = "id_sample") %>%
  group_by(id_sample) %>%
  summarise(x_bar = mean(price))

medias_amostra150 %>%
  ggplot(aes(x = x_bar)) +
  geom_histogram(
    aes(y = after_stat(density)), 
    binwidth = 2000, 
    fill = "orange", 
    col = "black"
    ) + 
  labs(
    title = "Distribuição das médias amostrais dos preços",
    subtitle = "2000 amostras com n = 150 cada"
  ) + 
  theme(text = element_text(size = 18)) + 
  stat_function(
    fun = dnorm, 
    args = list(mean = mean(medias_amostra150$x_bar), sd = sd(medias_amostra150$x_bar)),
    linewidth = 1, 
    color = "red"
    )

Como podemos ver no histograma acima, a forma da distribuição das médias para amostras de tamanho \(n=150\) é também normal, como esperado. Além disso, vemos também que nossa melhor estimativa para a média populacional continua sendo em torno de 180.000 dólares.

Comparando este caso com o da distribuição das médias dos preços para amostras de tamanho \(n = 15\), vemos que a maior diferença está na dispersão em torno da média populacional, que no presente caso é menor.

4. Das distribuições amostrais calculadas nos exercícios 2 e 3, qual tem menor dispersão? Se estamos interessados em estimativas que estão mais próximas do valor verdadeiro, preferiríamos uma distribuição com uma dispersão pequena ou grande?

Conforme já discutido no item acima, a distribuição amostral com menor dispersão é aquela com maior \(n\), que nesse caso é aquela com \(n = 150\). Se estamos interessados em fazer um inferência a respeito do valor da média populacional, por exemplo, deveríamos preferir amostras com maior \(n\), a fim de ter maior precisão.

medias_amostra_price = rbind(medias_amostra15, medias_amostra150) %>%
  select(-id_sample) %>%
  mutate(tamanho_amostra = factor(rep(c(15, 150), each = 2000)))

medias_amostra_price %>%
  ggplot(aes(x = x_bar)) + 
  geom_histogram(
    aes(y = after_stat(density)), 
    fill = "green", 
    col = "black", 
    binwidth = 2000
    ) +
  facet_wrap(~tamanho_amostra, ncol = 1, strip.position = "right") +
  theme(text= element_text(size = 18))

Conteúdos abordados neste laboratório

• Gráfico de pontos
• Correlação Linear
• Método de Mínimos Quadrados Ordinários
• Modelo de Regressão Linear Simples

Os dados

Vamos carregar os dados da temporada de 2011[^datasets]:

load(file = "mlb11.Rdata")

Além dos pontos marcados (runs), este conjunto de dados contém sete variáveis tradicionalmente utilizadas: vez ao taco (at_bats), rebatidas (hits), rebatidas para fora do campo (homeruns), média de rebatidas (bat_avg), eliminações (strikeouts), bases roubadas (stolen_bases), e vitórias (wins). Também foram incluídas três novas variáveis: percentual de alcance de base (new_onbase), percentual de potência (new_slug), e alcance de base mais potência (new_obs). Para a primeira parte da análise consideraremos as sete variáveis tradicionais. Ao final do laboratório, você trabalhará com as novas variáveis por conta própria.

Que tipo de gráfico você utilizaria para mostrar a relação entre runs e alguma outra variável numérica? Crie um gráfico dessa relação utilizando a variável at_bats como preditora. A relação parece ser linear? Se você soubesse o valor de at_bats de um time, você se sentiria confiante para utilizar um modelo linear para predizer o número de runs?

Para ter uma ideia da relação entre as duas variáveis, podemos usar um gráfico de pontos ou dispersão:

  ggplot(data = mlb11, aes(x = at_bats, y = runs)) + 
  geom_point(size = 2, color = "red") + 
  labs(
    title = "gráfico de dispersão",
    subtitle = "Pontos marcados como função de vezes ao taco"
  ) + 
  theme(text = element_text(size = 18))

Apesar da expressiva dispersão dos dados, é possível perceber, a partir do gráfico de pontos mostrado, que existe uma tendência de crescimento linear direto do número de pontos marcados (runs) com o número de vezes ao taco (at_bats). Portanto, dado um modelo linear devidamente ajustado aos dados por meio de algum algoritmo adequado de otimização (como mínimos quadrados ordinários, por exemplo), poderíamos sim utilizar esse modelo linear para predizer o número de pontos marcados por um time na temporada, usando a variável at_bats como preditora.

Se a relação parece ser linear, podemos quantificar a força da relação utilizando o coeficiente de correlação.

mlb11 %>%
  summarise(cor(runs, at_bats))

##   cor(runs, at_bats)
## 1           0.610627

Soma dos quadrados dos resíduos

Nesta seção você irá utilizar uma função interativa, plot_ss, para investigar o que queremos dizer com “soma dos quadrados dos resíduos”. Você deverá rodar esta função em seu console, não no documento markdown. Para rodar esta função é necessário que o conjunto de dados mlb11 esteja carregado em seu ambiente de trabalho.

Assim como utilizamos a média e o desvio padrão para resumir características importantes de uma única variável, podemos resumir a relação entre essas duas variáveis por meio de uma reta que melhor descreve sua associação. Utilize a seguinte função interativa para selecionar a linha que você acha que cruza a nuvem de pontos da melhor maneira.

plot_ss(x = mlb11$at_bats, y = mlb11$runs, showSquares = T)

Depois de executar esse comando, você será solicitado a clicar em dois pontos situados no gráfico para definir uma reta. Depois que você fizer isso, a reta que você especificou será mostrada na cor preta e os resíduos na cor azul. Perceba que há 30 resíduos, um para cada uma das 30 observações. Lembre-se que os resíduos são a diferença entre os valores observados e os valores preditos pela linha:

\[ e_i = y_i - \hat{y_i} \] A maneira mais comum de se fazer uma regressão linear é selecionar a linha que minimiza a soma dos quadrados dos resíduos. Para visualizar o quadrado dos resíduos, você pode rodar novamente o comando de geração do gráfico e adicionar o argumento showSquares = TRUE.

plot_ss(x = mlb11$at_bats, y = mlb11$runs, showSquares = T)

Note que o resultado da função plot_ss fornece a inclinação (coeficiente angular) e o intercepto da sua reta, bem como a soma dos quadrados.

O modelo linear

É bastante cansativo tentar obter a linha dos mínimos quadrados, ou seja, a reta que minimiza a soma dos quadrados dos resíduos, por meio de tentativa e erro. Ao invés disso, podemos utilizar a função lm no R para ajustar o modelo linear (também conhecido como linha de regressão).

m1 = lm(formula = runs ~ at_bats, data = mlb11)

O primeiro argumento da função lm é a fórmula, descrita como y ~ x. O que a função faz é obter o modelo linear de runs (pontos marcados) em função de at_bats (vezes ao taco). Aqui, cabe um destaque para o operador ~ (tilde), que é usado para separar as variáveis dependentes e independentes em um modelo. Neste caso, a fórmula diz: runs depende de at_bats. O segundo argumento especifica que o R deve buscar no banco de dados mlb11 as variáveis runs e at_bats.

A saída da função lm é um objeto que contém todas as informações que precisamos sobre o modelo linear que acabamos de ajustar. Podemos acessar essa informação utilizando a função summary.

summary(m1)

## 
## Call:
## lm(formula = runs ~ at_bats, data = mlb11)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -125.58  -47.05  -16.59   54.40  176.87 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2789.2429   853.6957  -3.267 0.002871 ** 
## at_bats         0.6305     0.1545   4.080 0.000339 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 66.47 on 28 degrees of freedom
## Multiple R-squared:  0.3729, Adjusted R-squared:  0.3505 
## F-statistic: 16.65 on 1 and 28 DF,  p-value: 0.0003388

Vamos analisar o resultado passo a passo. Primeiramente, a fórmula utilizada para descrever o modelo é apresentada no começo. Depois da fórmula você verá o resumo de cinco números dos resíduos. A tabela “Coefficients” (coeficientes) apresentada em seguida é fundamental; sua primeira coluna apresenta o intercepto-y do modelo linear e o coeficiente da variável at_bats. Com essa tabela, podemos descrever a reta de regressão de mínimos quadrados para o modelo linear:

\[ \hat{y_i} = -2789.2429 + 0.6305 \times at\_bats \] Uma última informação que abordaremos sobre a saída da função summary é o “Multiple R-Squared” (R-quadrado Múltiplo), ou de maneira abreviada, \(R^2\). O valor do \(R^2\) representa a proporção de variabilidade na variável resposta que é explicada pela variável explicativa. Neste modelo, 37,3% da variabilidade dos pontos runs é explicada pelas vezes ao taco at_bats.

Ajuste um novo modelo que utilize a variável homeruns para predizer runs. Utilizando as estimativas dos resultados do R, escreva a equação da reta de regressão. O que a inclinação (coeficiente angular) nos diz sobre a relação entre o sucesso de um time e seus home runs?

m2 = lm(formula = runs ~ homeruns, data = mlb11)
summary(m2)

## 
## Call:
## lm(formula = runs ~ homeruns, data = mlb11)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -91.615 -33.410   3.231  24.292 104.631 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 415.2389    41.6779   9.963 1.04e-10 ***
## homeruns      1.8345     0.2677   6.854 1.90e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 51.29 on 28 degrees of freedom
## Multiple R-squared:  0.6266, Adjusted R-squared:  0.6132 
## F-statistic: 46.98 on 1 and 28 DF,  p-value: 1.9e-07

Nesse caso, ajustando um modelo linear que utiliza a variável homeruns para predizer a variável runs, obtemos a seguinte equação para a reta de regressão:

\[ \hat{y_i} = 415.2389 + 1.8345 \times homeruns \] Ainda, o coeficiente angular nos diz que o número de pontos marcados (runs) cresce com o número de rebatidas fora do campo (homeruns) a uma taxa de 1,8345 vezes.

Para visualizarmos melhor a relação entre as duas variáveis, podemos plotar um gráfico de pontos:

ggplot(data = mlb11, aes(x = homeruns, y = runs)) +
  geom_point(size = 2, col = "red") + 
  labs(
    title = "Gráfico de dispersão",
    subtitle = "pontos marcados como função das rebatidas fora do campo"
  ) + 
  theme(text = element_text(size = 18))

Além disso, podemos ver também o coeficiente de correlação entre as variáveis runs e homeruns:

mlb11 %>%
  summarise(coef_cor = cor(runs, homeruns))

##    coef_cor
## 1 0.7915577

Predição e erro de predição

Vamos criar um gráfico de dispersão com a reta de mínimos quadrados para m1 disposta junto aos pontos.

ggplot(data = mlb11, aes(x = at_bats, y = runs)) +
  geom_point(size = 2, col = "red") +
  stat_smooth(method = "lm", formula = y ~ x, se = F) + 
  labs(
    title = "Modelo linear ajustado",
    subtitle = "vezes ao taco como preditora dos pontos"
  ) + 
  theme(text = element_text(size = 18))

Aqui nós estamos literalmente adicionando uma camada ao topo de nosso gráfico. A função stat_smooth cria uma reta ao ajustar um modelo linear. Ela também pode nos mostrar o erro padrão se associado à nossa reta, mas iremos omitir essa saída por enquanto.

Esta linha pode ser utilizada para predizer \(y\) em função do valor \(x\). Quando as previsões associadas aos valores de \(x\) são feitas além dos valores observados, nos referimos a elas como extrapolações e não é bem recomendado. Portanto, previsões realizadas à partir dos dados observados são mais confiáveis. Ela também é utilizada para computar os resíduos.

Se o gerente de um time visse a reta de regressão de mínimos quadrados e não os dados reais, quantos pontos (runs) ele prediria para um time com 5.579 vezes ao taco (at_bats)? Esse valor superestima ou subestima o valor real, e por quanto? Em outras palavras, qual é o resíduo para esta predição?

Utilizando os valores dos parâmetros ajustados de acordo com a amostra, \(\beta_0 = −2789.2429\) e \(\beta_1 = 0.6305\), além da equação do nosso modelo de regressão ajustado, podemos escrever:

p_marcados = as.double(subset(mlb11, at_bats == 5579, select = "runs"))

p_previstos = -2789.2429 + 0.6305*5579

cat(
  "Pontos previstos: ",
  p_previstos,
  "\nPontos marcados: ",
  p_marcados,
  "\nResíduo: ", 
  p_marcados - p_previstos, 
  "\nErro relativo absoluto: ", 
  abs((p_marcados - p_previstos)/p_marcados),
  "\n"
  )

## Pontos previstos:  728.3166 
## Pontos marcados:  713 
## Resíduo:  -15.3166 
## Erro relativo absoluto:  0.02148191

Escrita alternativa para o trecho acima:

p_marcados = mlb11 %>%
  filter(at_bats == 5579) %>%
  select(runs) %>%
  as.double()

p_previstos = -2789.2429 + 0.6305*5579

cat(
  "Pontos previstos: ",
  p_previstos,
  "\nPontos marcados: ",
  p_marcados,
  "\nResíduo: ", 
  p_marcados - p_previstos, 
  "\nErro relativo absoluto: ", 
  abs((p_marcados - p_previstos)/p_marcados),
  "\n"
  )

Conforme podemos ver na saída da linha de código acima, o número de pontos previsto pelo modelo em função do número de vezes ao taco, para o caso em que at_bats = 5579, é igual a 728, ao passo que o número de pontos observados para esse número de vezes ao taco é 713. Vemos, portanto, que o modelo superestima um pouco o valor realmente observado. Vemos também que a diferença relativa entre os dois valores é de cerca de 2%.

Diagnóstico do Modelo

Para avaliar se um modelo linear é confiável, precisamos verificar (1) a linearidade, (2) resíduos normalmente distribuídos, e (3) variância constante.

Linearidade: Você já verificou se a relação entre pontos (runs) e vezes ao taco (at_bats) é linear utilizando o gráfico de dispersão. Deveríamos também verificar essa condição utilizando um gráfico dos resíduos em função dos valores ajustados.

ggplot(data = m1, aes(x = .fitted, y = .resid)) +
  geom_point(col = "red", size = 2) +
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 1, col = "blue") + 
  labs(
    x = "valores ajustados",
    y = "Resíduos",
    title = "Modelo linear para runs x at_bats",
    subtitle = "Avaliação de viés nos resíduos"
  ) + 
  theme(text = element_text(size = 18))

Neste último gráfico, esperamos uma nuvem de pontos uniformemente distribuídos em torno da reta \(x = 0\), o que indicaria resíduos com média aproximadamente igual a zero e ausência de vieses.

Note que o nosso objeto m1 também funciona como um conjunto de dados, pois nele estão armazenados os valores ajustados (\(\hat{y}\)) e os resíduos. Observe tembém que nós utilizamos um código bem elegante aqui: depois de criar o gráfico de dispersão na primeira camada (primeira linha de código), nós adicionamos uma linha tracejada em \(y=0\) (para nos ajudar a checar se os resíduos estão distribuiídos em torno de 0), e também nomeamos os eixos de modo a tornar o gráfico mais informativo.

Há algum padrão aparente no gráfico de resíduos? O que isso indica sobre a linearidade da relação entre pontos (runs) e vezes ao taco (at_bats)?

Conforme podemos ver no gráfico de resíduos plotado, não há um padrão aparente (viés) na dispersão dos dados. Aparentemente, os resíduos estão dispostos aleatoriamente em torno de zero, o que é uma suposição incial em modelos de regressão. Contudo, a baixa quantidade de pontos dificulta um pouco a análise, já que alguns dos pontos estão muito afastados da reta tracejada.

Resíduos aproximadamente normais: Para verificar essa condição, podemos conferir o histograma dos resíduos:

ggplot(data = m1, aes(x = .resid)) +
  geom_histogram(
    aes(y = after_stat(density)), 
    binwidth = 20, 
    col = "black", 
    fill = "orange"
    ) + 
  stat_function(
    fun = dnorm, 
    args = 
      list(mean = mean(m1$residuals), sd = sd(m1$residuals)),
    col = "red", 
    linewidth = 1
  ) + 
  labs(
    x = "Residuals",
    title = "Distribuição dos resíduos"
  ) + 
  theme(text = element_text(size = 20))

ou um QxQ-plot (quantis dos resíduos vs quantis teóricos da normal padrão). Nesse caso, esperamos que se os quantis seguirem a mesma distribuição (normal) eles devem estar posicionados sobre uma reta com orientação positiva.

ggplot(data = m1, aes(sample = .resid)) + 
  geom_qq(col = "red", size = 2) + 
  geom_qq_line(linetype = "dashed", color = "blue", linewidth = 1) +
  labs(
    x = "Quantis teóricos",
    y = "Quantis amostrais",
    title = "Estudo da normalidade dos resíduos",
    subtitle = "q-q plot para os resíduos"
  ) + 
  theme(text = element_text(size = 20))

Note que a sintaxe para a construção deste gráfico é um pouco diferente do que você está aconstumado: nós indicamos sample para os resíduos aos invés de x (pois os dados provém de uma amostra, ou sample em inglês), e indicamos também o método estatístico geom_qq(), que significa “quantil-quantil” (QxQ), outro nome bem comum para este tipo de gráfico.

Conforme podemos obsersar no qXq plot e no histograma esboçado, os resíduos estão razoavlemente afastados da normalidade.

Variância constante: Essa condição pode ser verificada através da curva dos resíduos padronizados em função dos valores ajustados. Aqui, é importante ressaltar que essa condição também é conhecida como homocedasticidade.

Para verificar se um modelo de regressão linear ajustado atende aos pressupostos, podemos aplicar a função plot() à saída de lm(). Para entendermos qual será o resultado, podemos rodar o comando ?plot.lm.

• A função par() é usada para definir parâmetros de gráficos

• O argumento mfrow divide o gráfico conforme uma matriz \(n_r \times c_c\), plotando primeiro os quadros horizontais.

• A função plot(), quando aplicada a um objeto do tipo lm(), dispõe de 6 saídas:

1 - gráfico de dispersão para resíduos contra valores ajustados

2 - gráfico de localização da escala da raíz dos resíduos versus valores ajustados

3 - q-q plot para a distribuição normal

4 - gráfico de distâncias de Cook versus rótulos de linha

5 - resíduos versus alavancagem

6 - distâncias de Cook versus alavancagem

Por padrão, a saída apresenta as 3 primeiras e a quinta opções. Porém, outras podem ser selecionadas com o argumento which.

par(mfrow = c(2, 2))
plot(m1, which = c(1:4), pch = 16)

Olhando para o terceiro quadro acima, vemos que a hipótese de homocedasticidade aproximada dos resíduos está verificada.

Por sua conta

1. Escolha outra variável tradicional contida no banco de dados mlb11, diferente da variável at_bats, que você acha que poderia ser um bom preditor da variável runs. Crie um gráfico de dispersão das duas variáveis e ajuste um modelo linear. Visualmente, parece haver uma relação linear?

Vamos escolher a variável hits (rebatidas) como preditora. Deste modo, temos:

ggplot(data = mlb11, aes(x = hits, y = runs)) + 
  geom_point(col = "red", size = 2) + 
  stat_smooth(method = "lm", formula = y ~ x, se = T) + 
  labs(
    title = "Ajuste para pontos em função da média de rebatidas",
    subtitle = "modelo de regressão linear, ajuste linear"
  ) + 
  theme(text = element_text(size = 18))

mlb11 %>%
  summarise(cor(hits, runs))

##   cor(hits, runs)
## 1       0.8012108

Como podemos ver no resultado acima, existe uma aparente relação linear entre as variáveis runs e hits, que é corroborada pelo valor de \(\rho = 0.8012\).

2. Como essa relação se compara àquela existente entre runs e at_bats? Utilize os valores \(R^2\) do sumário dos dois modelos para compará-los. A variável que você escolheu parece predizer runs melhor do que at_bats? Como você justificaria sua resposta?

Para responder a essas perguntas, primeiro vamos obter os parâmetros da regressão:

m3 = lm(data = mlb11, formula = runs ~ hits)
summary(m3)

## 
## Call:
## lm(formula = runs ~ hits, data = mlb11)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -103.718  -27.179   -5.233   19.322  140.693 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -375.5600   151.1806  -2.484   0.0192 *  
## hits           0.7589     0.1071   7.085 1.04e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 50.23 on 28 degrees of freedom
## Multiple R-squared:  0.6419, Adjusted R-squared:  0.6292 
## F-statistic:  50.2 on 1 and 28 DF,  p-value: 1.043e-07

Conforme podemos concluir a partir dos valores do coeficiente de correlação entre as variáveis hitse runs, assim como a partir do valor do \(R^2\) múltiplo, a relação linear entre as duas variáveis se mostra mais forte do que aquela existente entre at_bats e runs. De fato, no atual caso em análise, o valor do \(R^2\) múltiplo é de aproximadamente 0.64, ao passo que no caso em que usamos a variável at_bats como preditora para runs, esse valor é de apenas 0.37. Portanto, concluímos que a variável hits (número de rebatidas) parece predizer melhor o número de pontos marcados pelo time na temporada do que o número de vezes ao taco (at_bats).

3. Agora que você pode resumir a relação linear entre duas variáveis, investigue a relação entre runs e cada uma das outras cinco variáveis tradicionalmente utilizadas no beisebol. Qual variável prediz melhor o valor de runs? Justifique sua conclusão utilizando métodos gráficos e numéricos já discutidos (para ser conciso, inclua apenas os resultados da melhor variável, não de todas as cinco).

Como já testamos 3 variáveis preditoras até agora (at_bats, home_runs e hits), vamos analisar as outras três e inserir cada resultado numa tabela:

m4 = lm(data = mlb11, formula = runs ~ bat_avg)
m5 = lm(data = mlb11, formula = runs ~ strikeouts)
m6 = lm(data = mlb11, formula = runs ~ stolen_bases)
m7 = lm(data = mlb11, formula = runs ~ wins)
m8 = lm(data = mlb11, formula = runs ~ new_onbase)
m9 = lm(data = mlb11, formula = runs ~ new_slug)
m10 = lm(data = mlb11, formula = runs ~ new_obs)

resultado1 = data.frame(
  var_predit = c("at_bats", "hits", "home_runs", "bat_avg", "strikeouts", "stolen_bases", "wins"), 
  r_square_mult = c(summary(m1)$r.squared, summary(m2)$r.squared, summary(m3)$r.squared, summary(m4)$r.squared, summary(m5)$r.squared, summary(m6)$r.squared, summary(m7)$r.squared), 
  cor_coef = c(cor(mlb11$runs, mlb11$at_bats), cor(mlb11$runs, mlb11$hits), cor(mlb11$runs, mlb11$homeruns), cor(mlb11$runs, mlb11$bat_avg), cor(mlb11$runs, mlb11$strikeouts), cor(mlb11$runs, mlb11$stolen_bases), cor(mlb11$runs, mlb11$wins))
  )

  kable(
    resultado1, 
    col.names = c("Variável Preditora", "R^2-múltiplo", "Coef. de Correlação")
    ) %>% 
    kable_styling(
      full_width = F, 
      bootstrap_options = c("striped", "hover", "condensed", "responsive")
      )

De acordo com a tabela apresentada acima, a melhor variável preditora para a variável runs, dentre as 7 listadas na tabela, é a variável bat_avg (média de rebatidas). Podemos ainda ver como fica a linha de regressão para esse caso:

ggplot(data = mlb11, aes(x = bat_avg, y = runs)) +
  geom_point(col = "red", size = 2) + 
  stat_smooth(method = "lm", formula = y ~ x, se = T) + 
  labs(
    title = "Ajuste para runs como função de bat_avg",
    subtitle = "modelo de regressão linear, ajuste linear"
  ) + 
  theme(text = element_text(size = 18))

4. Agora examine as três variáveis mais recentes. Essas são as estatísticas utilizadas pelo autor do filme Moneyball para predizer o sucesso de um time. De modo geral, elas são mais ou menos eficazes para predizer os pontos (runs) se comparadas às variáveis mais tradicionais? Explique utilizando evidências gráficas e numéricas apropriadas. De todas as dez variáveis que nós analisamos, qual parece ser a melhor preditora da variável runs? Utilizando as informações limitadas (ou não tão limitadas) que você conhece sobre estas estatísticas do beisebol, seu resultado faz sentido?

resultado2 = data.frame(
  var_predit = c("new_onbase", "new_slug", "new_obs"), 
  r_square_mult = c(summary(m8)$r.squared, summary(m9)$r.squared, summary(m10)$r.squared), 
  cor_coef = c(cor(mlb11$runs, mlb11$new_onbase), cor(mlb11$runs, mlb11$new_slug), cor(mlb11$runs, mlb11$new_obs)))

  kable(
    resultado2, 
    col.names = c("Variável Preditora", "R^2-múltiplo", "Coef. de Correlação")
    ) %>% 
    kable_styling(
      full_width = F, 
      bootstrap_options = c("striped", "hover", "condensed", "responsive")
      )

Conforme podemos ver na tabela acima, as estatísticas usadas no Filme Moneyball são melhores que todas as outras 7 tradicionalmente utilizadas para prever o número de pontos marcados por um time naquela temporada analisada. Dentre as novas estatísticas introduzidas, a variável new_obsé a que tem maior sucesso em predizer o número de pontos marcados por um time, considerando o \(R^2\) múltiplo e o coeficiente de correlação entre a variável runs e a respectiva variável preditora.

De fato, podemos ver como fica a linha de regressão para esse caso:

ggplot(data = mlb11, aes(x = new_obs, y = runs)) + 
  geom_point(col = "red", size = 2) +
  stat_smooth(method = "lm", formula = y ~ x, se = T) +
  labs(
    title = "Ajuste para runs como função de new_obs",
    subtitle = "modelo de regressão linear, ajuste linear"
  ) + 
  theme(text = element_text(size = 18))

5. Faça uma análise de diagnóstico do modelo de regressão a partir da variável que você escolheu como melhor preditora dos pontos (variável runs).

Para subsidiar a nossa análise, vamos aplicar a função plot() ao nosso modelo ajustado, m10, a fim de obter a avaliação dos resíduos do modelo.

par(mfrow = c(2, 2))
plot(m10, which = c(1:4), pch = 16)

Linearidade

Como podemos ver no gráfico acima, os resíduos não apresentam um viés claro aparente. Porém, a hipótese de que a média dos resíduos é zero parece estar violada.

Resíduos aproximadamente normais

Conforme podemos ver no qq-plot mostrado acima, os resíduos amostrais se ajustam razoavelmente bem a uma distribuição normal.

Variância constante

Com base no resultado do terceiro quadro mostrado acima, vemos que a hipótese de homocedasticidade é aproximadamente válida, sendo que a variância parece diminuir levemente com o aumento dos valores ajustados.

Finalmente, dados o gráfico do modelo ajustado e a análise de resíduos feita acima, podemos concluir que o modelo linear adotado (percentual de alcance de base como preditora para os pontos marcados) parece razoavelmente adequado. Contudo, devemos ter um pouco de cautela, já que duas das hipóteses utilizadas para a construção do modelo parecem ser ligeiramente violadas. Tal fato pode decorrer do reduzido número de pontos para análise (30).

Conteúdos abordados neste laboratório

• Teorema do Limite Central
• Distribuição Normal
• Distribuição normal padrão
• Intervalo de confiança

Especificando um semente: Nós iremos tomar algumas amostras aleatórias e calcular intervalos de confiança com base nessas amostras para este laboratório. Por esta razão, você deve especificar uma semente no topo do seu relatório. Caso tenha dúvidas sobre este tópico, retorne aos laboratórios anteriores.

set.seed(123)

Os dados

No quarto laboratório, “Distribuições Amostrais”, nós exploramos os dados populacionais das casas da cidade de Ames, Iowa. Vamos começar carregando este conjunto de dados.

load("ames.rda")

Neste laboratório, começaremos com uma amostra aleatória simples de 60 elementos da população. Perceba que o conjunto de dados contém informações sobre várias variáveis relativas às casas, mas para a primeira parte do laboratório focaremos no tamanho da casa, representada pela variável area.

n = 60
set.seed(123)
amostra = slice_sample(ames, n = 60, replace = T) 

cat(
  "Média Populacional:",
  mean(ames$area),
  "\nMédia Amostral:",
  mean(amostra$area)
)

## Média Populacional: 1499.69 
## Média Amostral: 1461.433

amostra %>%
  ggplot(aes(x = area)) + 
  geom_histogram(binwidth = 200, col = "black", fill = "orange") + 
  labs(
    y = "N",
    title = "Histograma para as áreas das casas em Ames",
    subtitle = "Amostra aleatória simples com 60 elementos"
  ) + 
  theme(text = element_text(size = 18))

Conforme podemos ver na saída dos comandos acima, a média populacional das áreas dos imóveis na cidade de Ames e a média da amostra aleatória simples de 60 casas que retiramos do conjunto original são muito parecidas. De fato, elas diferem por apenas 2,55%.

Esse é um resultado importante da estatística, que é uma consequência do teorema do limite central (TLC), o qual diz que as médias de amostras aleatórias simples de tamanho \(n\) terão distribuição aproximadamente normal quando \(n \rightarrow \infty\), e a média dessa distribição será igual à média da população subjacente. Por esse motivo, espera-se que outras pessoas que façam o mesmo procedimento amostral acima obtenham uma média amostral semelhante à média da população. Contudo, não é esperado que o valor que outra pessoa obtenha seja igual ao meu (se ela usar uma semente diferente durante o procedimento amostral), já que a média amostral é uma variável aleatória, que deve variar de uma amostra para outra.

Intervalos de Confiança

Retornemos, por um momento, à pergunta que motivou este laboratório: baseado nesta amostra, o que podemos inferir sobre a população? Baseado apenas nesta única amostra, a melhor estimativa da área habitacional das casas vendidas em Ames seria a média amostral, geralmente representada como \(\bar{x}\) (aqui denominaremos x_bar). A média amostral serve como uma boa estimativa pontual, mas seria interessante também comunicar quanta incerteza temos à respeito desta estimativa. Isso pode ser feito pelo uso de um intervalo de confiança.

Conforme pode ser acompanhado em qualquer livro-texto de estatística básica (MORETTIN; BUSSAB, 2017), um intervalo de confiança para a média de uma população com variância conhecida tem a seguinte forma

\[ \left( \bar{x} - z_{1-\alpha/2} \dfrac{\sigma}{\sqrt{n}}, \bar{x} + z_{1-\alpha/2} \dfrac{\sigma}{\sqrt{n}} \right). \]

Por questões didáticas, queremos trabalhar primeiramente com a distribuição normal. Contudo, é preciso deixar claro que a ideia de intervalo de confiança vale para qualquer distribuição de probabilidades. Por exemplo, se \(T\) for um estimador do parâmetro \(\theta\), e, se conhecemos a distribuição de \(T\), sempre será possível achar dois valores \(t_1\) e \(t_2\) que satisfaçam a relação (MORETTIN; BUSSAB, 2017)

\[ P(t_1 < \theta < t_2) = \gamma \],

sendo \(\gamma\) uma probabilidade, ou seja, \(0 < \gamma < 1\). No caso específico das distribuições normais e da t-Student, os dois quantis são simétricos, já que a prória distribuição é simétrica, em ambos os casos.

Uma vez que queremos usar as medidas amostrais e não calculamos a verdadeira variância da população, em nosso exemplo iremos supor normalidade dos dados e considerar que o tamanho da amostra é suficientemente grande para utilizarmos o desvio padrão amostral \(s\) no lugar de \(\sigma\). Assim não teremos que, obrigatoriamente, utilizar a distribuição t-Student. Nosso intervalo de confiança terá, então, a forma

\[ \left( \bar{x} - z_{1-\alpha/2} \dfrac{s}{\sqrt{n}}, \bar{x} + z_{1-\alpha/2} \dfrac{s}{\sqrt{n}} \right). \]

Aqui, é preciso fazer um parêntesis: Teoricamente, ao usar o desvio padrão amostral \(S\) em vez do desvio padrão populacional \(\sigma\), que é o que fazemos na prática, pois não conhecemos os parâmetros populacionais quando estamos fazendo inferência, a distribuição da estatística do teste muda para uma distribuição t de Student. Porém, também sabemos que, quanto maior o número de graus de liberdade da t de Student (que está relacionado ao tamanho da amostra), mais essa distribuição se aproxima de uma normal padrão. Por isso, esta aproximação está sendo usada aqui, já que o autor não queria introduzir a distribuição t de Student neste momento do estudo.

Neste ponto do curso você já deve se sentir confortável para calcular a média e o desvio padrão amostrais no R, e sabemos que o tamanho da amostra é 60. Desse modo, o único ponto que nos falta é determinar o valor crítico (quantil) associado a um nível de confiança dado. Lembre-se que níveis de confiança e quantis não são equivalentes. Por exemplo, um nível de confiança de 95% se refere a 95% da região central da distribuição, e o valor crítico associado a essa área na distribuição normal é o quantil de 97,5% da distribuição. Como a distribuição é simétrica, o outro quantil é o negativo deste.

Intervalo de confiança de 95% para a média, incluindo os respectivos quantis da normal padrão.

Nós podemos encontrar o valor crítico para um intervalo com confiança de 95% usando

z_gamma = qnorm(0.975)

que é o valor 1,96 com o qual você já está bastante familiarizado.

Vamos finalmente calcular o intervalo de confiança:

amostra %>%
  summarise(inf = mean(area) - z_gamma*sd(area)/sqrt(n), sup = mean(area) + z_gamma*sd(area)/sqrt(n))

## # A tibble: 1 × 2
##     inf   sup
##   <dbl> <dbl>
## 1 1348. 1575.

Ou, como é mostrado nos livros-textos:

cat("IC(mu, 95%) = ",
    "[",
    mean(amostra$area) - z_gamma*sd(amostra$area)/sqrt(n),
    ",",
    mean(amostra$area) + z_gamma*sd(amostra$area)/sqrt(n),
    "]"
    )

## IC(mu, 95%) =  [ 1348.355 , 1574.512 ]

Recapitulando: mesmo que não saibamos como a população inteira se distribui, temos 95% de confiança de que a verdadeira média do tamanho das casas em Ames se encontre entre os valores inf (limite inferior do intervalo de confiança) e sup (limite superior do intervalo de confiança). Contudo, existem algumas condições que precisam ser atendidas para esse intervalo ser válido.

Para que esse resultado seja válido, a população (que nesse caso é composta pelas áreas das casas na cidade de Ames) precisa ter distribuição normal, já que, para determinar um intervalo de confiança para a média com variância desconhecida, precisamos considerar a variável aleatória \(\frac{\sqrt{n}(\bar{X}-\mu)}{S}\), a qual terá distribuição t de Student se a população já tiver distribuição normal. Além disso, estamos considerando o fato de que, quando \(n \rightarrow \infty\), a distribuição t de Student com \(n\) graus de liberdade tende para uma normal padrão. Nesse caso, com \(n = 60\), a aproximação é razoável, dando uma diferença de cerca de 2% no quantil relativamente ao resultado exato.

Níveis de Confiança

O que significa 95% de confiança? Em poucas palavras, significa que, se repetirmos esse procedimento e determinarmos a média amostral um número grande de vezes, cerca de 95% deles conterá o valor da média da população.

Em nosso caso, temos o luxo de sabermos a verdadeira média populacional, uma vez que temos os dados da população inteira. Vamos calcular este valor para determinar se nosso intervalo de confiança realmente captura a média. Vamos armazenar este valor com o nome param (abreviação para parâmetro populacional), e nomeá-lo mu.

param = ames %>%
  summarise(mu = mean(area))
param

## # A tibble: 1 × 1
##      mu
##   <dbl>
## 1 1500.

Assim, vemos que o nosso intervalo de confiança de 95% nesse caso contêm a média \(\mu\) da população. Além disso, se muitas pessoas fizessem esse mesmo procedimento amostral feito aqui, espera-se que cerca de 95% dos intervalos de confiança gerados contenham o valor real de \(\mu\).

Utilizando o R, vamos criar várias amostras para aprender um pouco mais a respeito de como as médias amostrais e os intervalos de confiança variam de uma amostra para outra.

Eis o esboço do processo:

1. Obter uma amostra aleatória.

2. Calcular a média e o desvio padrão da amostra.

3. Utilizar estas estatísticas para calcular um intervalo de confiança.

4. Repetir estas etapas 50 vezes.

Nós podemos alcançar este objetivo por meio da função rdply do pacote plyr. As seguintes linhas de comando tomam 50 amostras aleatórias simples de tamanho n da população (lembre-se que definimos n=60 mais cedo), e calcula os intervalos de confiança para cada uma das amostras.

set.seed(123)
ic = rdply(50, slice_sample(ames, n = n, replace = T), .id = "Amostra") %>%
  group_by(Amostra) %>%
  summarise(inf = mean(area) - z_gamma*sd(area)/sqrt(n), sup = mean(area) + z_gamma*sd(area)/sqrt(n))

Vamos dar uma olhada nos cinco primeiros intervalos:

ic %>%
  slice(1:5) %>%
  kable() %>%
  kable_styling(
    full_width = F, 
    bootstrap_options = c("striped", "hover", "condensed",
                          "responsive")
    )

Em seguida iremos criar um gráfico bastante ilustrativo: representaremos cada intervalo de confiança que captura a média pela cor azul e cada intervalo que não captura a média pela cor vermelha. Note que capturar a média significa que o limite inferior (inf) está abaixo do valor da média e o limite superior (sup) está acima do valor da média. O primeiro passo na construção deste gráfico é atualizar a variável ic, adicionando a seguinte informação: ele captura ou não a média? Faremos isso por meio da função mutate.

ic = ic %>%
  mutate(captura_mu = ifelse(inf < param$mu & sup > param$mu, "sim", "não"))

A função ifelse é nova. Ela carrega três argumentos: primeiro uma declaração lógica, em segundo lugar vem o valor que desejamos apresentar se a resposta à declaração lógica for verdadeira, e em terceiro lugar vem o valor que desejamos apresentar se a resposta à declaração por falsa.

Agora possuímos toda a informação necessária para a criação do gráfico, mas necessitamos reorganizar nossos dados um pouquinho para facilitar o “desenho” do gráfico. Especificamente, nós queremos organizar os dados em um novo banco de dados onde cada linha represente um limite relacionado a um intervalo, ou seja, uma variável chamada ic_id irá identificar a qual intervalo nos referimos (isso tanto para o inf como para o sup), a variável ic_limites irá guardar os limites inferior e superior de cada intervaldo indexado pela ic_id e a variável captura_mu irá nos dizer se aquela linha corresponde a um limite que captura ou não a verdadeira média populacional. Podemos alcançar esse objetivo por meio do código abaixo:

ic_dados = data.frame(
  ic_id = c(1:50, 1:50), 
  ic_limites = c(ic$inf, ic$sup), 
  captura_mu = c(ic$captura_mu, ic$captura_mu)
  )

E, finalmente, podemos criar o gráfico usando o seguinte:

ic_dados %>%
  ggplot(aes(x = ic_limites, y = ic_id, group = ic_id, color = captura_mu)) + 
  geom_point(size = 3) +
  geom_line() +
  geom_vline(xintercept = param$mu, color = "darkgray") +
  labs(
    x = "IC(95%)",
    y = "Amostra",
    title = "Intervalos de confiança para a média",
    subtitle = "amostras de tamanho n = 60"
  ) + 
  theme(text = element_text(size = 18))

Qual a proporção de intervalos de confiança que capturam a verdadeira média da população? Essa proporção é exatamente igual ao nível de confiança? Se não, explique o motivo. Não se esqueça de incluir seu gráfico à sua resposta.

Conforme podemos ver no gráfico acima, apenas 1/50, ou seja, 2% dos intervalos de confiança gerados, não contêm o verdadeiro valor da média. Teoricamente, esse valor deveria ser igual a 5%, já que estamos usando um intervalo de confiança de 95%. Contudo, o processo de amostragem possui natureza aleatória e, portanto, a depender do número de amostras que estivermos considerando, o número de intervalos que não conterão a verdadeira média naturalmente irá oscilar, mas certamente sua fração correspondente não se distanciará muito de \(1 - \gamma\), que nesse caso é igual a 5%.

Por sua conta

1. Escolha um nível de confiança de sua preferência, desde que não seja 95%. Qual é o valor crítico apropriado?

Vamos escolher, por exemplo, dois níveis de confiança tipicamente presente nos livros-textos ou nos trabalhos científicos, que são os de 90% e 99%:

qnorm(c(0.95, 0.995))

## [1] 1.644854 2.575829

Nesse caso, geramos os valores críticos (quantis), já bastante conhecidos na literatura, que são \(z_1 \approx 1.64\) e \(z_2 \approx 2.57\).

2. Calcule 50 intervalos de confiança utilizando o nível de confiança que você escolheu na questão anterior, e construa um gráfico similar ao que acabamos de criar, e calcule a proporção de intervalos que incluem a verdadeira média da população. Como se compara essa proporção com o nível de confiança que você escolheu para os intervalos? Não se esqueça de incluir seu gráfico à sua resposta.

set.seed(123)
z_gamma2 = qnorm(0.95) # Quantil de referência para o intervalo de 90%
ic2 = rdply(50, slice_sample(ames, n = n, replace = T), .id = "Amostra") %>%
  group_by(Amostra) %>%
  summarise(
    inf = mean(area) - z_gamma2*sd(area)/sqrt(n), 
    sup = mean(area) + z_gamma2*sd(area)/sqrt(n)) %>% 
  mutate(captura_mu = ifelse(inf < param$mu & sup > param$mu, "sim", "não"))

ic_dados2 = data.frame(
  ic_id = c(1:50, 1:50), 
  ic_limites = c(ic2$inf, ic2$sup), 
  captura_mu = c(ic2$captura_mu, ic2$captura_mu)
  )

ic_dados2 %>%
  ggplot(aes(x = ic_limites, y = ic_id, group = ic_id, color = captura_mu)) +
  geom_point(size = 3) +
  geom_line() +
  geom_vline(xintercept = param$mu, color = "darkgray") +
  labs(
    x = "IC(90%)",
    y = "Amostra",
    title = "Intervalos de confiança para a média",
    subtitle = "amostras de tamanho n = 60"
  ) + 
  theme(text = element_text(size = 18))

Conforme podemos ver no gráfico gerado, apenas três intervalos em cinquenta, ou seja, 6% dos intervalos de confiança gerados não incluem a média da população. Novamente, se repetirmos esse experimento com um grande número de amostras, esperamos que aproximadamente 10% dos intervalos de confiança gerados não incluam a verdadeira média da população, já que o nível de confiança foi fixado em 90%. Portanto, o nível de significância, \(1-\gamma = 10\%\), representa o percentual de amostras que não incluirão a média da população.

Conteúdos abordados neste laboratório

• Intervalo de confiança
• Teste t para a média de uma população supostamente normal
• Teste t pareado e não pareado para a diferença entre duas médias

Nascimentos na Carolina do Norte

Em 2004, o estado da Carolina do Norte, Estados Unidos, disponibilizou um grande conjunto de dados contendo informações sobre os nascimentos registrados no estado. Esse conjunto de dados é útil para pesquisadores que estudam a relação entre hábitos e práticas de gestantes e o nascimento de seus filhos. Nós trabalharemos com uma amostra aleatória das observações deste conjunto de dados.

Análise Exploratória

Carregue o conjunto de dados nc em seu espaço de trabalho (OPENINTRO, [s. d.]).

load("nc.RData")

Temos dados de 13 variáveis diferentes, algumas qualitativas e outras quantitativas. O significado de cada variável é o seguinte:

Primeiro, vamos fazer um breve resumo desse conjunto de dados:

cat("Número de Registros:",
    nrow(nc),
    "\nResumo Estatístico:\n"
    )

## Número de Registros: 1000 
## Resumo Estatístico:

summary(nc)

##       fage            mage            mature        weeks             premie   
##  Min.   :14.00   Min.   :13   mature mom :133   Min.   :20.00   full term:846  
##  1st Qu.:25.00   1st Qu.:22   younger mom:867   1st Qu.:37.00   premie   :152  
##  Median :30.00   Median :27                     Median :39.00   NA's     :  2  
##  Mean   :30.26   Mean   :27                     Mean   :38.33                  
##  3rd Qu.:35.00   3rd Qu.:32                     3rd Qu.:40.00                  
##  Max.   :55.00   Max.   :50                     Max.   :45.00                  
##  NA's   :171                                    NA's   :2                      
##      visits            marital        gained          weight      
##  Min.   : 0.0   married    :386   Min.   : 0.00   Min.   : 1.000  
##  1st Qu.:10.0   not married:613   1st Qu.:20.00   1st Qu.: 6.380  
##  Median :12.0   NA's       :  1   Median :30.00   Median : 7.310  
##  Mean   :12.1                     Mean   :30.33   Mean   : 7.101  
##  3rd Qu.:15.0                     3rd Qu.:38.00   3rd Qu.: 8.060  
##  Max.   :30.0                     Max.   :85.00   Max.   :11.750  
##  NA's   :9                        NA's   :27                      
##  lowbirthweight    gender          habit          whitemom  
##  low    :111    female:503   nonsmoker:873   not white:284  
##  not low:889    male  :497   smoker   :126   white    :714  
##                              NA's     :  1   NA's     :  2  
##                                                             
##                                                             
##                                                             
## 

Com a ajuda da função glimpse(), vamos ver quais dados são numéricos e quais são do tipo categórico:

glimpse(nc)

## Rows: 1,000
## Columns: 13
## $ fage           <int> NA, NA, 19, 21, NA, NA, 18, 17, NA, 20, 30, NA, NA, NA,…
## $ mage           <int> 13, 14, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 16,…
## $ mature         <fct> younger mom, younger mom, younger mom, younger mom, you…
## $ weeks          <int> 39, 42, 37, 41, 39, 38, 37, 35, 38, 37, 45, 42, 40, 38,…
## $ premie         <fct> full term, full term, full term, full term, full term, …
## $ visits         <int> 10, 15, 11, 6, 9, 19, 12, 5, 9, 13, 9, 8, 4, 12, 15, 7,…
## $ marital        <fct> married, married, married, married, married, married, m…
## $ gained         <int> 38, 20, 38, 34, 27, 22, 76, 15, NA, 52, 28, 34, 12, 30,…
## $ weight         <dbl> 7.63, 7.88, 6.63, 8.00, 6.38, 5.38, 8.44, 4.69, 8.81, 6…
## $ lowbirthweight <fct> not low, not low, not low, not low, not low, low, not l…
## $ gender         <fct> male, male, female, male, female, male, male, male, mal…
## $ habit          <fct> nonsmoker, nonsmoker, nonsmoker, nonsmoker, nonsmoker, …
## $ whitemom       <fct> not white, not white, white, white, not white, not whit…

Vamos agora investigar a idade típica das mães e pais na Carolina do Norte. Comecemos construindo histogramas, boxplots e verificando o resumo de cinco números para cada caso:

histogram = nc %>%
  ggplot(aes(x = mage)) + 
  geom_histogram(
    col = "black", 
    fill = "pink", 
    binwidth = 1) + 
  labs(
    title = "Mother's age",
    x = "age",
    y = "N"
  ) + 
  theme(text = element_text(size = 18))

box_plot = nc %>%
  ggplot(aes(y = mage, x = "")) + 
  geom_boxplot(
    fill = "pink",
    outlier.colour = "red",
    alpha = 0.95,
    width = 0.4
  ) + 
  labs(
    x = "",
    y = ""
  ) +
  theme(text = element_text(size = 18))

cowplot::plot_grid(histogram, box_plot)

histogram = nc %>%
  ggplot(aes(x = fage)) + 
  geom_histogram(col = "black", fill = "skyblue", binwidth = 1) + 
  labs(
    title = "Father's age",
    x = "age",
    y = "N"
  ) + 
  theme(text = element_text(size = 18))

box_plot = nc %>%
  ggplot(aes(y = fage, x = "")) + 
  geom_boxplot(
    fill = "skyblue",
    outlier.colour = "red",
    alpha = 0.95,
    width = 0.4
  ) + 
  labs(
    x = "",
    y = ""
  ) +
  theme(text = element_text(size = 18))

cowplot::plot_grid(histogram, box_plot)

cat("Resumo de 5 números para a idade da mãe:\n")

## Resumo de 5 números para a idade da mãe:

summary(nc$mage)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      13      22      27      27      32      50

cat("Resumo de 5 números para a idade do pai:\n")

## Resumo de 5 números para a idade do pai:

summary(nc$fage)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   14.00   25.00   30.00   30.26   35.00   55.00     171

cat("\nVariabilidade da idade da mãe:\n")

## 
## Variabilidade da idade da mãe:

nc %>% 
  summarise(dq = IQR(mage, na.rm = T), s = sd(mage, na.rm = T))

##   dq        s
## 1 10 6.213583

cat("\nVariabilidade da idade do pai:\n")

## 
## Variabilidade da idade do pai:

nc %>% 
  summarise(dq = IQR(fage, na.rm = T), s = sd(fage, na.rm = T))

##   dq        s
## 1 10 6.763766

Note que sd(fage) e IQR(fage) não retorna resultados válidos. O resumo de 5 números indica que existem 171 nascimentos para os quais a idade do pai não foi informada. Por padrão, a maioria das funções do R não retorna resultados válidos quando o dado está ausente. Isto pode ser resolvido adicionando-se o argumento na.rm = TRUE na chamada da função.

Introdução à inferência

Suponha que queremos testar a hipótese de que a idade média para as mães que dão à luz na Carolina do Norte é 26,5 anos. Calculando “manualmente”, temos(MORETTIN; BUSSAB, 2017):

# Dados amostrais

xbar = mean(nc$mage)
std = sd(nc$mage)
n = length(nc$mage)

# Teste de hipótese -- Calculando o p-valor

dof = n - 1 #Graus de liberdade
std_err = std/sqrt(n) #erro padrão do estimador
test_stat = (xbar - 26.5)/std_err #Estatística do teste
p_value = 2*pt(-abs(test_stat), df = dof)
p_value

## [1] 0.01108826

#Construindo um intervalo de confiança de 95% para a média:

q = qt(0.975, df = dof)
IC = xbar + c(-1, 1)*q*std_err
IC

## [1] 26.61442 27.38558

Como o nosso p-valor é menor que o nível de significância \(\alpha\) (0,05), temos evidências suficientes para rejeitar a hipótese de que a idade média das mulheres que dão à luz na Carolina do Norte é de 26,5 anos. Observe que 26,5 anos não está contido no intervalo de confiança de 95% construído. Os testes de hipóteses para a média com nível de significância \(\alpha\) são logicamente equivalentes a um intervalo de confiança de \(100(1 - \alpha)\%\) - essa é a grande sacada.

Vamos agora fazer o mesmo procedimento com o uso da função t.test():

# Mother's mean age -- HT and 95% CI

t.test(
  x = nc$mage, 
  alternative = "two.sided", 
  mu = 26.5, 
  conf.level = 0.95
  )

## 
##  One Sample t-test
## 
## data:  nc$mage
## t = 2.5446, df = 999, p-value = 0.01109
## alternative hypothesis: true mean is not equal to 26.5
## 95 percent confidence interval:
##  26.61442 27.38558
## sample estimates:
## mean of x 
##        27

Suponha agora que queiramos testar se a idade média dos pais da NC é 30 anos. Use \(\alpha = 0.01\). Construa também um intervalo de confiança de 99% para a média. Calcule “manualmente” e usando t.test(). Dica: quando calculamos “manualmente”, os valores ausentes podem ser um problema, de modo que primeiro precisamos extrair e armazenar as informações úteis. A função t.test(), por outro lado, já cuida disso automaticamente.

# Dados Amostrais

# f_age = nc %>%
#   select(fage) %>%
#   drop_na()

xbar = mean(nc$fage, na.rm = T)
std = sd(nc$fage, na.rm = T)
n = sum(!is.na(nc$fage))

# Teste de hipótese - calculando o p-valor

dof = n - 1 #Graus de liberdade
std_err = std/sqrt(n) #Erro padrão do estimador
test_stat = (xbar - 30)/std_err #Estatística do teste

2*pt(-abs(test_stat), df = dof)

## [1] 0.2766452

# Construindo um intervalo de confiança de 99% para a média:

q = qt(0.995, df = dof)
IC = xbar + c(-1, 1)*q*std_err
IC

## [1] 29.64923 30.86223

Fazendo o mesmo cálculo usando a função t.test(), temos:

t.test(
  x = nc$fage, 
  mu = 30, 
  alternative = "two.sided",
  conf.level = 0.99
  )

## 
##  One Sample t-test
## 
## data:  nc$fage
## t = 1.0886, df = 828, p-value = 0.2766
## alternative hypothesis: true mean is not equal to 30
## 99 percent confidence interval:
##  29.64923 30.86223
## sample estimates:
## mean of x 
##  30.25573

Suponha que um pesquisador quer testar a hipótese de que a média da idade dos pais na NC é diferente da média da idade das mães que deram à luz - suponha um \(\alpha\) de 0,01. Este é um teste para a diferença entre as duas médias populacionais, o que nos obriga a perguntar se os dados para os dois grupos (mães e pais) são pareados ou não. Claramente, os dados são pareados, uma vez que as idades do pai e da mãe guardam relação de interdependência. Fazendo o cálculo primeiro “manualmente”, temos:

# Teste pareado - apenas é necessária a diferença nas idades

age_diff = nc$mage - nc$fage
summary(age_diff)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## -23.000  -5.000  -2.000  -2.653   0.000  12.000     171

# Dados amostrais

xbar = mean(age_diff, na.rm = T)
std  = sd(age_diff, na.rm = T)
n = sum(!is.na(age_diff))

# Teste de hipótese - calculando o p-valor

dof = n - 1 #graus de liberdade
std_err = std/sqrt(n) #erro padrão do estimador
test_stat = (xbar - 0)/std_err #estatística do teste
2*pt(-abs(test_stat), df = dof)

## [1] 1.504608e-59

# Construindo um intervalo de confiança de 99% para a média

q = qt(0.995, df = dof)
IC = xbar + c(-1, 1)*q*std_err
IC

## [1] -3.040107 -2.265080

E usando a função t.test(), nosso resultado sai muito mais rápido:

# Opção 1: Calcular a diferença entre as duas variáveis externamente e
# inserí-la na função t.test()

t.test(x = age_diff, alternative = "two.sided", mu = 0, conf.level = 0.99)

## 
##  One Sample t-test
## 
## data:  age_diff
## t = -17.673, df = 828, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 99 percent confidence interval:
##  -3.040107 -2.265080
## sample estimates:
## mean of x 
## -2.652593

# Opção 2: Faça a própria função calcular as diferenças, utilizando
# o argumento paired = TRUE

t.test(
  x = nc$mage, 
  y = nc$fage, 
  paired = T, 
  alternative = "two.sided",
  mu = 0,
  conf.level = 0.99
  )